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ABSTRACT 

The normal mode expansion technique is applied to the transformed 
monoenergetic integral transport equation to develop a solution for the 
rotationally invariant and axially infinite, critical, two-region cylin- 
der with a finite outer reflector boundary. The model assumes isotropic 
scattering and identical neutron mean free paths in the core and reflec- 
tor regions. The solution in terms of singular integral equations is 
obtained by applying a completeness theorem found for the singular 
eigenfunctions. Numerical results for a variety of core and reflector 
multiplying properties and reflector thicknesses are presented and com- 
pared with the results of other methods. The completeness inherent in 
this solution and the high precision in the numerical calculations pro- 
vide results which may be used as analytic standards for this problem. 

An example of this type of application is given in a study of approxi- 
mations inherent in the neutronic design analysis of a small, fast- 
neutron-spectrum reactor concept proposed as a space power source. 

Using the highly precise critical dimensions for the case which most 
closely approximates this reactor, investigations were made of the 
reactivity effects of angular quadrature type and order and two- 
dimensional geometrical models used in the discrete -ordinates transport 
analysis of this concept. 
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CHAPTER I 
INTRODUCTION 


1.1 Purpose and Scope 

The purpose of this dissertation is to present in detail the 
development of a highly precise transport solution for the radially 
reflected critical cylinder and to demonstrate how the results from 
this solution can be used as analytic standards in evaluating approxi- 
mations inherent in numerical transport treatments employed in reactor 
design analysis. 

A complete solution for the rotationally invariant and axially 
infinite two-region critical cylinder with a finite outer reflector 
boundary is obtained by applying the singular eigenfunction expansion 
technique to the transformed monoenergetic integral transport equation. 
The model assumes isotropic scattering and identical neutron mean free 
paths in the core and reflector regions. Numerical results for a 
variety of core and reflector multiplying properties and reflector 
thicknesses are presented and compared with the results of other methods. 
The critical dimensions and the neutron density distribution for one 
of these cases are then used as analytic standards in evaluating dis- 
crete angular segmentation transport programs used in reactor design. 

Three aspects of the numerical programs were studied. 

1. Angular quadrature order (number of segmentations in the 
angular variable) 

. Type of angular quadrature (direction cosines and weights 
chosen by various prescriptions) 


2 
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3. The step-boundary approximation inherent in X-Y geometrical 
representations of circular boundaries. 

1.2 Background and Dissertation Organization 

Neutron transport as a function of position, energy, angle and 
time is generally assumed to be described by the Boltzmann equation. 

In Chapter II this equation is reduced to the time -independent, one- 
dimensional, monoenergetic form. In 1960, Case 1 obtained a complete 
solution to the reduced Boltzmann equation in one -dimensional plane 
geometry in terms of singular eigenfunctions. This solution is briefly 
described in Chapter II. Concluding Chapter II is a demonstration of 
the equivalence between the integral form of the transport equation 

and the Boltzmann equation. 

2 

In 1963 Mitsis obtained exact solutions for the critical sphere 
and for the critical infinite cylinder by transforming the monoenergetic 
integral transport equation and applying the singular eigenfunction 

3 

expansion method of Case. Case and Zweifel have shown a more general 
treatment of the same problems by demonstrating a replication property 
of the kernel of the integral transport equation. The replication 
method has been extended by Gibbs^ to obtain solutions in arbitrary 
convex geometry. For the critical sphere and cylinder, Gibbs special- 
ized his general solution to duplicate the results of Mitsis. 

5 

Lathrop and Leonard have suggested that numerical results from 
the Mitsis solution for the critical cylinder could be used to investi- 
gate the accuracy of two-dimensional discrete -ordinates angular quadra- 
ture sets. However, a recent bibliography 6 of neutral particle transport 
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theory does not contain any reference to such results and initial 
investigation under the present study was directed towards determining 
the efficacy of the Mitsis solution for the hare critical cylinder. 

This investigation found that the solution as formulated by Mitsis 
is not convergent. However, a converging solution could be obtained 
by reformulating a function used in separating the outer boundary 
condition into singular and nonsingular parts and by redefining the 
continuum expansion coefficient. Another variation of the solution 
for the bare critical cylinder presented in Chapter III from that of 
Mitsis is a more straightforward, though equivalent, treatment of the 
singular eigenfunctions. The advantages of this new singular eigen- 
function treatment are discussed below. 

The development of the solution for the radially reflected critical 

cylinder is presented in Chapter IV. The solution is based upon the 

integral transform approach developed for the bare core solution by 

Mitsis. The same approach has been taken by Smith and Siewert^ and 
8 

by Leuthauser in developing solutions for the reflected critical 
sphere. Both of these solutions assume identical neutron mean free 
paths in the core and reflector regions. Smith 9 has demonstrated the 
complicated form of the transformed integral equation when this assump- 
tion is not made. To reduce the complexity of the problem, the same 
assumption is imposed on this solution for the reflected critical 
cylinder. However, the multiplying properties of the two media are 
allowed to differ. In the context of the monoenergetic model, the 
identical mean free path assumption is not a severe restriction. 
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especially for fast energy spectrum reactors where total neutron cross 
sections are of the same order of magnitude. 

Transformation of the integral transport equation in two-region 
cylindrical geometry results in integral equations for the neutron 
density in each region. Contained in these expressions are kernels 
made up of stuns of integrals of modified Bessel functions over the 
spatial variable. These kernels are separated out and defined as 
pseudo neutron distribution functions. The pseudo neutron distribution 
functions are shown to satisfy the same integro-differential equation 
as arises in the bare core case with the same centerline and outer 
boundary conditions. Additional boundary conditions are continuity 
of the pseudo neutron distribution function and its spatial derivative 
at the core -reflector interface. 

Solutions in terms of modified Bessel functions and singular 
eigenfunctions (called pseudo eigenfunctions because they are functions 
of the transformed variable) are found for the integro-differential 

o 

equation by the separation of variables technique. When Mitsis 
introduced these pseudo eigenfunctions, he initially considered the 
full range on the eigenvalues. Then he observed that the pseudo 
eigenfunctions correspond to the sum of the Case plane geometry eigen- 
functions and showed completeness by the theorem proven by Case. ^ 

In the present solutions for both the bare and reflected systems, 
advantage is taken of the evenness of the pseudo eigenfunctions and 
of the dispersion function for the eigenvalues to consider only those 
eigenvalues in the positive half-range. A full-range completeness 
theorem for the pseudo eigenfunctions is proven in Appendix B and its 
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easy extension to the half -range is shown. This procedure provides 
a basis for developing the solution entirely in tewns of the eigen- 
functions found, without using a decomposition into the Case plane 
geometry eigenfunctions as done by Mitsis in obtaining the bare core 
solution. A particular advantage of the present approach is the lack 

of dependence on the half-range plane geometry eigenfunction complete- 

3 

ness theory which requires the calculation of X-functions. Also 
inherent in the proof of completeness are several very useful ortho- 
gonality and normalization relations for the pseudo eigenfunctions. 
These relations are developed in Appendix A. 

Substitution of the pseudo eigenfunctions into the boundary 
conditions and application of the orthogonality-normalization relations 
results in two coupled iterative sequences with which the expansion 
coefficients can be calculated to any desired degree of accuracy. The 
critical condition arises from the reduction of the outer boundary 
condition to a Fredholm integral equation for one of the reflector 
expansion coefficients. It corresponds directly to an auxiliary con- 
dition required in the proof of completeness. 

Chapter V contains a description of the numerical techniques used 
to obtain the various functions and parameters appearing in the solu- 
tions developed in Chapters III and IV. The accuracy of these tech- 
niques is drawn from comparison with tabulated values. Then follows 
a demonstration of how the precision of the results varies with the 
order of numerical quadrature used in evaluating the integral terms 
appearing in the solutions. Finally the numerical results for a wide 
variety of cases are presented and compared with the results of other 
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methods. The bare core results are compared with the results of other 
analytic solutions found in the literature. The reflected core results 
are compared with the results of high order numerical calculations. 

Chapter VI contains a demonstration of how the results from the 
highly precise solution for the reflected cylinder can be used as 
analytic standards, A case is chosen which, in the context of the 
model, most closely approximates the Advanced Power Reactor^ concept 
being studied at the Lewis Research Center. Then, using the exact 
critical dimensions, the discrete angular segmentation programs used 
in the design of this concept are employed to study the three effects 
listed under the Purpose and Scope of this study (1. 1) . 

Chapter VII contains the conclusions and recommendations drawn 
from this study. The conclusions pertain both to the analytic solu- 
tions found and to the accuracies of the various numerical design 
approximations studied. Based on these conclusions, recommendations 
are made as to order of approximation required for desired design 
accuracies and as to where indicated sources of error might be reduced. 



CHAPTER II 


MONOENERGETIC TRANSPORT THEORY 


2.1 The Boltzmann Equation 

The assumptions under which the Boltzmann equation describes 
neutron transport in a medium free of independent sources are listed. 

1. The neutron acts as a point particle which travels in a 
straight line with constant speed between neutron-nuclei 
interactions. 

2. The probability of neutron-neutron interactions is much less 
than the probability of neutron-nuclei interactions and is, 
therefore, ignored. 

3. The neutron population is sufficiently large so that statisti- 
cal fluctuations may be ignored. 

4. Secondary neutrons are produced at the time and position of 
the primary neutrons, 

5. The total cross section, a(r,E), describing the probability 
of a neutron -nucleus interaction per unit path length, is a 
function of energy and position only. 

With these assumptions, the Boltzmann equation for a medium free 
of independent sources is written as'' - '*' 


1 6 |_ 
v 6 1 


(r,ft,E,t) + 0_*V§ (r,Q,E,t) + a (r,E) $ (r,ft,E,t) 



(r,E ’ ) f(r;0',E'-»n,E) f (r^^E^tjdfi 'dE' 


( 2 . 1 ) 
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where the neutron angular flux, $(r,Q,E,t), Is a function of position, 
angle, energy and time; v is the neutron speed; and cr (r,E) x 
f (r 1 ,E ’->0, E) gives the total probability of neutron transfer from 
n',E' to Q,E. 

The present study is restricted to the greatly simplified, mono- 
energetic, stationary form of the Boltzmann equation. Elimination of 
the energy and time dependence in Eq. (2.1) yields 

( 2 . 2 ) 

where c(r) is defined as the mean number of neutron secondaries per 
collision and enters from the operation 

X f(r;n , ,E , -jfi,E)dE'= c(r) f“(r;fi'-0 
_ ' 


fi-V§(r,fi) + a (r) $ (r,n) = a (r)c(r) / f(r;ft'-£3) $(r,0')dD 

*Ti' 


where c(r), o (r) and ffr;^'-^) are one group, spectrum averaged values. 

In order to demonstrate Case’s singular eigenfunction technique, 
Eq. (2.2) is written for infinite plane geometry with dependence on 
one coordinate only. In planar geometry we have the relations 

5$ A 

= M- r — , where u, = Q*z, 

— 6 z ’ — * 

and also dQ ’ = 2ndp', -1 ^ u ^ 1, 


With these substitutions, Eq. (2,2) becomes 

.1 

&4 . .. . . . > r , , , . , v- / 

H- 


5 $ f 1 

(z,|l) + Cf (z)f (z,|i) = o(z)c(z)2n J f (z;fi ’-sfi) $ (z,h ' )dp. ' , 

-1 


Eq. (2.3) is further simplified by the restriction to isotropic 
scatterings in the laboratory system. In this instance. 


(2.3) 
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f(z }£'-«) = 1/4 TT. 


Finally, the spatial variation is expressed in neutron mean free paths 


through the definition: 


x = I o (z ' )dz ' 


from which follows ^ = a (z) ~ . 


With the above substitutions, Eq. (2.3) becomes 


h (x,|i) + #(x. 


,H) = t(x,l 


h ')du ' . 


Eq. (2.4) is the stationary, monoenergetic, one -dimensional, Boltzmann 
equation for a source-free, isotropic, homogeneous medium in plane 
geometry for which Case 1 developed the singular eigenfunction solution. 
2.2 Case's Singular Eigenfunction Method 

The solution to Eq. (2.4) is sought through the separation of 
variables technique. Thus we assume a general solution of the form 


$(x,m,) = X(x)<p(|i) 


Substitution of Eq. (2.5) into Eq. (2.4) yields 


h<p(h) + X(x)cp(n) = | X(x) I cp(|i, 1 )d|j,’ 


After collecting like variables in Eq. (2.6) we obtain 


1 dX(x) c_ 

X(x) dx ~ 2jicp( 


-1 ^ 
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Next, each side of Eq. ( 2 . 7 ) is set equal to the separation constant, 
— l/ u, to obtain the two equations: 


and 


dX(x) X(x) 0 
dx u 


H 


= f cp(^ 


)<V 


The solution of Eq. (2.8) is immediately written as 

X(x) = (constant )e~ X ^ v . 


To obtain the solution for Eq. (2-9), we first normalize 

cp(u)d|x = 1 

after which Eq. (2.9) becomes 

(u - u)cp(u) = 4r 

Eq. (2.11) yields a family of nonsingular and singular solutions 
depending on the value of u. For u/[-l,l], we have two discrete solu 
tions given by 

f ^ _ C U Q 
“ 2 u 0 Tp 

where ± u 0 are the roots of the dispersion function 

A(u) = 1 - cutanh _ ^(l/u) = 0 . 

For ue[-l,l], we have a continuum of solutions 

cp u (p.) = f p ” + X(u)&(u-m,) 

where \(u) is the function of Case defined by 

X(u) = 1 - cutanh 1 (u), 



( 2 . 8 ) 

(2-9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

( 2 .lU) 


(2.15) 
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and P indicates that the Cauchy principal value is taken for integrals 
involving cp (4). Collecting Eqs. (2.10, 2.12, 2.1b), we write the general 
solution of Eq. (2.^) as 


$(x,4) = a 0+ cp 0+ (4)e‘ X//u °+ a 0 _ep 0 _(4)e x / u °+ f ^ A(v)q> u (\i)e X / U du. 


(2.16) 


-1 S |l s 1 


Case 1 proved a completeness theorem for the above expansion which is 
constructive in the sense that it provides a method for obtaining the 
expansion coefficients. He also proved a partial range completeness 
theorem. Since these proofs are rather involved and are well documented, 
we do not include them here. However, an analogous proof for the eigen- 
functions developed in Chapter III is presented in Appendix B. 

2.3 The Integral Transport Equation 

Attempts by Mitsis to apply Case’s method for solving the Boltzmann 
equation in other than plane geometry were unsuccessful. However, suc- 
cessful applications were obtained in spherical and cylindrical geometry 
after making suitable transformations of the integral form of the trans- 
port equation. Here we wish to show the equivalence of the Boltzmann 
equation and the integral transport equation. First the monoenergetic, 
stationary Boltzmann equation as it appears in Eq. (2.2) is rewritten as 

n • v$(r ,Q) + a(r)«(r,n) = q(r,o) 
where the source term, q(r,n), is given by 

q(r,Q) = cr(r)c(r) / f(r;n' - f>)§ (r,0' )dfj’ . 


(2.17) 
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With the position vectors as shown below, 



we define r = £p + sn and write Eq. (2-17) as 


|| (£o + en,n) + a (rp + sn)$(£p + sq,q) « <i(£o + sn,n) • 

Operation on Eq. (2.18) with the integrating factor 

exp |^* ct(£p + s*Q)ds*j 


(2.18) 


results in 


^^expjjTo^ +s'n)ds , ]< 


r s 

*(£o +h 0»Q)| = cr(lo + s'n)ds J 

( 2 . 19 ) 

as can be shown by applying Leibnitz's Rule. In Eq. (2.19) we replace 
s by s' and s' by s", then operate by the integral over s' from -« 
to s obtaining 

s 

expf y* ^ (r o +s"n)ds'J$(r 0 +s'n,n) ■L expj^” o (r^+s’'n)ds'^Jq(r 0 +s 'Q_,Q)ds ' » 

“co W (2.20) 

In evaluating the left hand side of Eq. (2.20) at the lower limit we 
assume that 

lim $ (r 0 + 8*0,0) = 0. 


q(r 0 + s0,0) 


S 1 — ► —oo 


Transposing the exponential term on the left hand side of Eq. (2.20), 
dividing the range on the integral in the transposed exponential from 
-co to s 1 and from s’ to s, and, after cancellation, we obtain 



13 


§(r. 


+ s M fl)ds" 


+ s'n,n)ds' . (2.21) 


+ sn,n) = L exp “L/! c( -° 

Now recalling that r = Tp + so, r^ = £ - sQ, we write Eq. (2.21) as 

$(r,0) = J exp^-^ o(r + [s" - s]o)ds'^ q(r + [s’ - s]0,n)ds’ . (2.22) 

Next, inserting new variables p" = s - s" and p 1 = s - s’, Eq. (2.22) 
becomes 


*(r,Q) = I expT- £ a(r -p"n)dp"Jq(r -p’n,n)dp' . 


(2.25) 


Returning to the more conventional displacement variable, we let 
p* = s and p" = s’ in Eq. (2.25) obtaining 

CO S 

$(r,n) = I expj -j a(r - s’ 0 )ds , J q(r - sQ,n)ds. (2.24) 

Upon imposing isotropic scattering on the source term, thereby eliminating 
its dependence on O’ as it appears in Eq. (2.17), Eq. (2.24) becomes 

CO S 

$(r,Q) = */o expjj. / cr(r - s 'n)ds’J cr(r - sO)c(r - sq)$ (r - sO)ds. (2.25) 
In order to eliminate the angular dependence, we operate on Eq. (2.25) 


with the integral over Q. as follows: 

00 r s 

y* = exp yj' ~ s -) ds 


n o 


a (r - sO)c(r - sO)$(r - sQ_)dsdO, 

(2.26) 


let r' = r - sfi, and recognize that the volume element dsc£2. = 
dr'/|r-r'|^, and use the definition for optical length, T(|r'-r|) = 


r'-r 


J ^ ofr-s'njds', to obtain 

exp[ -t( |r ' - rj )] 


§ (£) = &J 


c (r')c(r ! )§ (r ' )dr 1 , 


(2.27) 


r - r ’ 
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After expressing the volume element dr' as d^r' and changing the 
spatial variable to mean free paths, Eq. (2.27) becomes 

*(r) = ijjff c(r’)$(r’)d 3 r T . ( 2 . 28 ) 

Finally we rewrite Eq. (2.28) in terras of the neutron density by 
multiplying through by the average speed, v, and using p(r) = v$(r), 
to obtain 

^TF 0( - )o( - )d3r '‘ (2 ' 29) 

Eq. (2.29) is the stationary, monoenergetic, integral transport equation 
for the neutron density in a source free, homogeneous medium with 


isotropic scattering. 



CHAPTER III 


SOLUTION FOR THE BARE CRITICAL CYLINDER 


3.1 Transformation of the Integral Equation in Cylindrical Geometry 
The monoenergetic integral transport equation, Eq. (2.29), is 


written for a single region as follows : 

, -|r-r' 

/ \ c / p(r')e 

»<i> — 

J r< I £ " £' I 2 


d 3 r ', 


where c, the mean number of neutron secondaries per collision, is a con- 
stant. We wish to apply this equation to the bare cylindrical geometry 
shown in Fig. 1. In Fig. 1 we have represented the position r* by the 
cylindrical coordinates (t, a, z) and have located r at (r, 0 , o) . We 
observe from Fig. 1 that 

|r - r' I 2 = x 2 + z 2 = r 2 + 1 2 - 2rtcos(a) + z^ . (3 


Expressing the volume integral over r’ in cylindrical coordinates. 


Eq. (3.1) becomes 


p(r) = ^ / tp(t)dt 


2n a 

f m f 

* / 0 -'-oo 


- „/x g + z s , 

X a + Z a 


where R represents the radius of the outer boundary and the limits on 
z correspond to the axially infinite cylinder. 


The integral over z is reduced by first defining I as 

Z 


co 

2 f exp [- x 

x s / 1 

•/o 


+ (z/x ) 2 


+ (z/x 


Then I can be written as 
z 


CO 


00 

/ exp [- u + (z/x) 2 ] du. 


(3-4b) 
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Substituting v= z/x, Eq. (3*4b) becomes 


CO QO 

;/ du / 


exp u 71 + 

«A + va 


Now by defining t' = ^/l + v 2 , we can express Eq. (3* 4c) as 

I= i/ au /£^L. (; 

Z X x *1 yt’ s -1 

The second integral in Eq. (3«4d) corresponds to the Laplace transform 
of F (t’) = (t' s -1)”^^, which is equal to K o (u), modified Bessel 
function of the second kind; therefore, 

CO 

I = - f K (u)du. (I 

z x J o v 

X 

Finally, we make the substitution, p, = x/u, into Eq. (3*4e) obtaining 

i 

1=2 f dp,. (I 

Z J 0 M- 2 

Substitution of Eq. (3*4f) into Eq. (3*3) and expressing x as in 


Eq. (3-2) results in 


2tt 1 


n ^ r / 

'<>■) - b J t / * / K c(- 

O o O' 


+ t 2 - 2rtcos (a) \ du 


The integral over a in Eq. (3-5) is performed by applying the addition 
theorem for the modified Bessel function of the second kind; namely, 


+ 1 ; 


■ 2rtcos "(a 


)\ Y ina) r S t 

’/ ^ * j K n (t/p,)l n (r/ w ), r g t 


and noting that 


e ina da = 2 Z> n = ° 


0, n ^ 0 « 



18 


The result of this procedure is an integral equation for the neutron 
density distribution tfith integrals over the radial variable, t, and 
the transformed variable, p,. 


pO) 


= c 


f ^\f K o (r/n)l o (t/|iHp(t)dt 


+ y K o (t/n)l o (r/(j,)tp(t)dtJ (3-7) 

3.2 The Pseudo Neutron Distribution Function 

2 

Adopting the nomenclature of Mitsis, we define the kernel appearing 
in Eq. (3*7) as a pseudo neutron distribution function, $(r,|i), related 
to the neutron density such that 


where 


_L 

p( r ) = f ^7$^ 

o ^ 

= cjy*K o (r/fi)l o (t/n)tp(t)dt 

R 

+ / K 0 {t/n)l o (r/n)tp(t)dt 


(3-8) 


(3-9) 


Next, we wish to show that $(r,|i) obeys an integrodifferential equation 
made up of first and second derivative terms. In talcing the first 
derivative, we use the relations: 


k;<z) - (z) = i^z). 

Use of Leibnitz's Rule and the above relations yields the first derivative 
of $(r,|j,) with respect to r as 

SSfr.n) c I / % K 1 (r/^)l o (t/^)tp(t)dt + (t/n)l 1 (r/^)tp(t)dt . (3.10a) 

6r m- \J r J 
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The second derivative is found by continuing the above procedure and 
using the additional relations: 

I{(z) = I D (z) - \ I 1 (z);K|(z) = -K o (z) - \ K x (z); 

12 

and the Wronskian for the modified Bessel functions: 

K c (z)l 1 (z) +K 1 (z)l o (z) = l/z. 


Use of the above relations yields the second derivative of $(r,p,) as 


' £ /[ 


K 0 (r/^) + £ K^r/fi) 


I 0 (t/n)tp(t)dt 


+ J ^ K CJ (t/jj,)^I o (r/|j 1 ) - ft I^(r/^,)^jtp(t)dt 


- cp(r). 

Now substituting from Eqs- (3-9)> (5- 10a) and (3-10b), we form the 
equation: 


+ 1 &*(r,n) 
&r r ftr 


= 


y|jc o (r/pi) +£ K 1 (r/n)jl o (t/n)tp(t)dt + 
j/ r K o (t/ p )j^ I o ( r /4) -ft ^(r/njjtpftjdt - 

K 1 (r/n)l 0 (t/ u )tp(t)dt + f K 0 (t/n)l 1 (r/|x)tp(t)dtj- 


c p (r) + 


Al 


K Q (r/ M.)l Q (t/ (j,)tp (t )dt + j£ K 0 (t/n)l 0 (r/n)tp(t)dt 


- 


Canceling like terms on the right hand side of Eq. (3* 10c) an d substi 
tuting Eq* (3*8) for p(r) we obtain 


(3.* 10b) 


(3*10c) 
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+ iSllM) _ 1 #( r ,n) = -c / (5-10d) 

6r 2 r br |x 3 U- 

Having shown the integro-differential equation that $(r,(j) obeys, now 
we need to establish the boundary conditions satisfied by $(r,p,). 

Inspection of $(r,p,) as defined in Eq. (5*9) at the core centerline, 
r = o, with the consideration that I (0) = 1 and that the limit of 


t K (t/n) = o as t -»o, establishes the boundary condition: 

$(o,p) finite. 


(3.11a) 


The outer boundary condition is constructed by first substituting r = R 
into Eq. (3.9) to determine $(R,|J>) as 

»(R,H) = c f K (R/n) 1 (t/^)tp(t)dt. (3. 

•'O 

Next, substituting r = R into Eq. (3«10a) yields for the first deriva- 
tive term 


n 

= - £ / K 1 (R/ f x)l o (t/^)tp(t)dt. 

r=R ^ Jo 


(3-llc) 


Inspection of Eqs. (3.11b and c) as combined below establishes the 
outer boundary condition as 


K C (R /n) 


« (r,^) = 0. 


(3. Hd) 


3«3 Singular Eigenfunction Expansion and Solution 

In developing the eigenfunctions, we assume a solution of Eq. (3*10d) 
of the form 


$ (r,^) = R (r )M (p.) 
v V V 


(3.12) 
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2 

Substitution of Eq. (3*12) into Eq. (3«10d) and using l/ v as the 
separation constant yields 

d 2 Ry(r) 1 dR v (r) 1 


dr^ 


+ r dr 


v 


R (r) = 0 


and 


/I 1\ /*>'> 

(i? * ^■) M v (|x) = ° J dp ’* 


The general solution of Eq. (3*13) is 


R v (r) = al c (r/v) +PK o (r/y). 


To find the eigenfunctions and eigenvalues of Eq. (3 -14), we first 
normalize by setting 

/y 


m> •) 

— <V = i 


after which we have 


(v a - jj, 3 )M^(|ju) = cvV, ° - M- - 1 


Here we find it convenient to make the definition 


T)>) = m v (m.)/ b 2 , 

where we call our pseudo eigenfunction. Substitution of Eq. 

into Eq* (3*17) yields one discrete solution as 

C 

= vS -Ja ’ ^0,1] 

where v 0 is obtained from the dispersion function, Eq. (2.13). For 
ve[0,l], a continuum of solutions is obtained as 

\(») = cP + X (v)b(v - |i) 


(3.13) 

(3-14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 
(3.18) 

(3.19) 


( 3 . 20 ) 
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where P denotes that the Cauchy principal value is taken for integrals 

involving T)^(|j,), and \(v) is the function of Case'*' given by Eq. (2.15). 

2 

When Mitsis introduced these pseudo eigenfunctions, he initially 
considered both positive and negative ranges on the eigenvalues, and 
then observed that the pseudo eigenfunctions are related to the plane 
geometry eigenfunctions of Case, Eqs. (2.12) and (2.l4), in the following 
manner: 

Vli) = ^0+^) +C P Q _(^) 

and 

T) v (m.) = «P V (|*> +cp_ v (m,)- 

Mitsis went on to complete his bare core solution in terms of the plane 
geometry eigenfunctions. Therefore, the Mitsis solution depends on 
half-range plane geometry eigenfunction completeness theory which 
involves the calculation of X- functions. Here we take advantage of 
the evenness of the pseudo eigenfunctions and the dispersion function 
to consider only those eigenvalues in the positive half range. In 
Appendix B we prove a full-range completeness theorem for the pseudo 
eigenfunctions and show its easy extension to the half range. This 
procedure provides a basis for developing the solution entirely in 
terms of the pseudo eigenfunctions. In addition to the theoretical 
consistency of the present approach, the resulting solution provides 
a simpler sequence for obtaining numerical results than the approach 
taken by Mitsis. The numerical results from both approaches are 
equivalent. Furthermore, the proof of completeness for the pseudo 
eigenfunctions involves orthogonality and normalization relations 
developed in Appendix A which are very useful in obtaining the solution 
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for the two-region problem. Efforts to obtain a two-region solution 
with the Mitsis approach, that is, using the plane geometry eigen- 
functions, were unsuccessful. 

Proceeding with the bare core solution, the general solution for 
§(r,|i) is obtained by substituting Eqs. (5-15)? (3*19) and (3-20) into 
Eq. (3*12) yielding 


0 s r £ R $(r,p,) = [a o I o (r/v 0 ) +b Q K o (r/v o ) ]^ 2 ll o ((i) 

0 s n ^ 1 

1 

+ /"[A(v)l (r/v) +B(v)K Q (r/v) ]|i 2 T] v (|i)dv. (3*21) 

*o 

For a critical system, c, the mean number of neutron secondaries 

per collision, must be greater than unity. In this instance, as 

given by Eq. (2.15) will be purely imaginary. Thus, in Eq. (5*21), 

r o( r /v o ) becomes J 0 ( r /|v Q |) and K Q (r/v 0 ) is equal to \ i(J 0 (r/ |v q |) + 

iY(r/jv o |))« With the singular nature of Y Q (r/|v o |) and K Q (r/v) at 

r = o, b and B(v) must be zero to satisfy the boundary condition, 
o 

Eq. (3«lla). Eq. (3.21) for the pseudo distribution function is now 
rewritten as 

1 

*(r,n) = s-qJ (r-/ | | ) M. 2 Tl o ( m. ) +/a(v)I (r/v)|i 2 Tl (|x)dv. (3-22) 

o v 

The next step in the solution is to substitute §(r,p,) as defined 
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2 

Multiplying Eq. (3.23) by R/p , we obtain 


Jk 1 (h/m.)J 0 (b/|v 0 |) k 0 (b/ |1 ) j 1 (ii/|v 0 |) 

I 5 1^1 

1 

+ J' A(v)q.(v»M>)‘ll (M-) d v = 0, 
o 


(3.24) 


where 


q(\un) = ^K 1 (r/m,)I 0 (r/v) +^C 0 (b/h)I 1 (b/v). 


At this point we wish to separate Eq* (5*24) into singular and 
nonsingular parts. This is effected by defining the function: 

1 


H(v^) - 


V - jJL 


i o (r/m-) | 

TTr 7^T <l( ''’ l * ) " 1 l 


(3.25a) 


or 


< — I, (R/v) | 

H(v.li) = ^ ^ K 0 (E/u)l 0 (E/n) - !(• 0 - 25 b} 


Considering the behavior of the modified Bessel functions as shown in Fig. 

2 and the asymptotic expressions, Eqs. (5.10) and (5.11), pg. 55, it can 

be seen that H(v,mO is a bounded function. This procedure closely parallels 

2 

the work of Mitsis. However, the function of Mitsis which corresponds to 

H(v,p,) is unbounded. The reformulation of H(v,p.) as presented here is sug- 

4 

gested in the general solution of Gibbs for arbitrary convex bodies. 

Since v and ^ range over the same values, it can be shown with 
the use of the Wronskian that H(v,p), as defined in Eq. (3.25b), is 
indeterminate. Use of L’ Hospital's rule yields H(v,v) as 

H(v,v) = K q (r/v) |i^(r/v)/I 0 (r/ v) -I q (r/v)| (3*25c) 


which is also a bounded function. 
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Next we rewrite Eq. (3 -25a) as 

i 0 (r/v) 

q.(wM-) = J {l + (v - ji)H(v,u)| 


(3.25d) 


and substitute for q(v>^) in Eq. (3*24) obtaining 

} l(R/v) . 

j A(v) T {l + (v - n)H(v,p,)} Tl v ((i)dv 

(k 1 (r/h)j o (r/|v 0 |) k 0 (r/ m )j 1 (r/|v 0 |) 


KT 


Ra oV^' 


(3.26a) 


Now returning to the pseudo eigenfunction as defined in Eq. (3.20) we 
note that 

11 1 
f (v-n)T) (fi)dv - f (v-n)TcP +\ (v)B(v-n)ldv = c f dy 

(3.26b) 

Using this observation, multiplying through by I q (r/^), and transposing 
terms, we complete the separation of Eq. (3* 26a) into singular and non- 
singular parts as 


-L 

/a’(v)1>)* = (m-) 


(3.27a) 


where 


§’(n) = I (r/h) 


k o (r/^)j 1 (r/|v o |) 


k 1 (r/ m .)J 0 (r/|v 0 |) 




O 'O 


and 


_ c f A T (v)H(v,ki)v s dv 

-4 v + n 


A’(v) = A(v)I q (r/v) . 


(3.27b) 

(3.27c) 
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The discrete expansion coefficient, a Q , is chosen to correspond to 
some arbitrary power level and is thereby set equal to unity. The remain 
ing problem is to find a solution for the continuum expansion coefficient 
A T (v). Eq. (3.27a) does not provide a closed form solution for A’(v) 
because of the nondegenerate Fredholm term appearing in the definition 
of $'(n). However, an iterative solution for A' (v) can be constructed, 
provided that a free expression for A T (v) can be found. Now Eq. (3.27a) 
is the equivalent half-range form of Eq. (B.2a) in the proof of complete- 
ness for the pseudo eigenfunctions. Therefore, the desired free expres- 
sion for A'(v) corresponds to Eq. (B.15) which we write as 

1 

A '(v) = N^JJ ' Ch)dh (3*28) 

where N(v) is the normalization function for the continuum eigenfunctions 
Eq. (A- 31). The criticality condition is the auxiliary condition Eq. 
(B.lla) in the completeness proof, expressed over the half range as 


Eqs. (3.27a, b, c), 
bare critical cylinder, 
the following scheme : 

*>) = I 0 (k/m,)RT| o ( |J u)J— - X l — | ( 3 • 30b ) 


_L 

f 1 * 2*1 

•l + 


(4) 


c / ' |j dj! = o . 

'o o I 


(3-29) 


(3.28), and (3.29) form the solution for the 
Successive approximations are constructed by 


a;(ul) = o 


(3. 30a) 


IK (R/u)J, (R/I vI) K,(b/u)J. (R/ lv. I) 
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1 


o 

(3.30c) 

f A A+i (v)H(v^)v 3 

«' , (u) = $'(») -C / -=£= — dv 

n+l vp/ o v J J v + Li 

0 

(3.30d) 

r *•;*<■*> , „ 

c 4 ' 

(3- 30e) 


The critical core radius is found by conducting a radius search to 
satisfy the criticality condition, Eq. (3-30e), for a convergent itera- 
tive sequence for A'(v), Eqs. (3*30a through 3-30d). After the critical 
radius and converged expansion coefficients are obtained, the neutron 
density distribution can be calculated as 

1 

p(r) = a o J o (r/|v 0 |) + J A' (v)dv, (5-31) 

o 

which we have obtained by operating on Eq. (3.22) as in Eq. (3.8), 
substituting A'(v) from Eq. (3.27c), and using the normalization, Eq. 
(3.16). Numerical results for the bare critical cylinder are presented 


in Chapter V. 



CHAPTER IV 


SOLUTION FOR THE REFLECTED CRITICAL CYLINDER 

4.1 Formulation of the Equations 

The solution for the two-region problem generally follows that of 
the one-region problem. However, here we have the added complexity of 
an intermediate boundary separating the core region C^ > 1, and the 
reflector region, Cg < 1. We begin by rewriting Eq. (2.29), the mono- 
energetic integral transport equation describing the neutron density 
in an isotropically scattering medium as 


p(r) 



(4.1) 


where c(r_' ) denotes the mean number of neutron secondaries per collision 
and distances are measured in mean free paths. By assuming a constant 
neutron mean free path throughout, we are able to apply this equation to 
the two region cylindrical geometry shown in Fig- 3- In Fig. 3 we have 
represented r' by the cylindrical coordinates (t, a, z) and have located 
r at (r, o, o). We observe from Fig. 3 that 

■ . ,2 2 2 2.2 _, . . 2 

|r - r 1 ] - x + z = r +t -2rtcos(a) +z 


Expressing the volume integral over r’ in the cylindrical coordinates, 
Eq. (4.1) becomes 


R/- 


2tt 


p(r) = T~ f f act J eXP ^ 8 z 5 +z£ J dz, (4.2) 




Figure 3. - Two region cylindrical 
geometry. 
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where c(t) = C^, 0 S t £ R^; c(t) - C 2 , R^ ^ t i Rg, and the limits on 
z correspond to. the axially infinite cylinder. 

The integral over z in Eq. (4.2) is reduced in the same manner 
as in Chapter III, Eqs. (3.4a through 3«4f). Under this procedure 
Eq. (4.2) becomes 

p(r) = £ / o(t)tp(t)dt f 'da f k (^ . (4.3) 

*o o o ' ^ ^ 


The integral over a in Eq. (4.3) is performed by applying the addition 
theorem for the modified Bessel function of the second kind, 



+ t s - 2rtcos(q) 
H 



K n (r/n)l n (t/^) 
^ n (t/ u)l n ( r / M- ) 


r S t 
r i t > 


and noting that 



2rr, 

0, 


n = 0 
n^O * 


The result of this procedure is an integral equation for the neutron 
density in each region given by 


l , T 

p l( r ) - / <C X /K 0 (r/n)l (t/n)tp(t)cLt 


olrgR, 


R 1 

\f K (t/n)l (r/n)tp(t)dt 

C 2 f K 0 (t/^)l o (r/^)tp(t)dtl , 

R 1 ‘ 


(4.4a) 
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P 2 (r) 




K 0 (r/n)l 0 (t/n)tp(t)dt 


\ S r 5 ^ 


C 2 y*K o (r/ 4 )l o (t/n)tp(t)dt + C 2 J* K o (t/n)l o (r/n)tp(t)dt|. 


(4.4b) 


4.2 The Pseudo Neutron Distribution Functions 

As in the bare core solution, we again adopt the nomenclature 
2 

of Mitsis by defining the kernels appearing in Eqs. (4.4a & b) as 
pseudo neutron distribution functions related to the neutron densities 
such that , 


where 


and 



f *1 (*> 4 ) 



Pl( r ) = / ^2 

(4.5a) 

and 

1 



, N / #a( r ^4) J 

p 2 r = J d *‘> 

(4.5b) 


r 

$ 1 (r,u) = C 1 / K (r/n)l (t/n)tp(t)dt 

•A 


o S r S R. 


R, 


1 

+ C3 f K 0 (t/ u )l 0 (r/n)tp(t)dt 

+ c 2 / K o (t/^)l o (r/ p )tp(t)dt, 
*1 

J 2 (r, u ) =C f K (r/j*)l (t/n)tp(t)dt 

‘'O 

R, IrlR. 

12 r 

+ C 2 ^*K o (r/ M ,)l o (t/ (1 )tp(t)dt 


(4.6a) 


h 

C 2 / K 0 ( t /|J-)l 0 (r/|J.)tp(t)dt. 

Jr 


(4.6b) 
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Next we wish to show that both $^(r,u) and ^(r^p,) °4ey the sajne 
integro-differential equation, Eq. (3-10d), obeyed by the pseudo neutron 
distribution function in the bare core solution. Toward this end we 
take the first derivative of $^(r,y,) and ^(r,^) by using the relations 
K'(z) = -K, (z) and I'(z) = I..(z), and Leibnitz's rule to obtain 

O 1 0-*- 


6$ 


and 


i— . = - T / Mr/iOiJt/iOtpWat 

jT 

+ ^ J" K o (t/ M ,)l 1 (r/n)tp(t)dt 

c r 2 

+ — J K o (t/n)l 1 {r/u)tp(t)dt, 


& r 


. a. Z" 1 K (r/„)I (t/„)tp(t)dt 
^ Jo 

r 

A(r/u)l(t/,)t P (t)dt 

* R x 

c r 2 

+ — f K (t/ 4 )l 1 (r/^)tp(t)dt. 
h J r 


(4.?a) 


(4.?b) 


We find the second derivatives by continuing the above procedure, 
using the additional relations 

I{(z) = I D (z) - \ \{*)> K{(z) = - K q (z) - | K^z) 
and applying the Wronskin for the modified Bessel functions 

K^zjl^z) + K^z)^) = l/z 

to obtain 
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e> 2 $i(r,ix) 


6 r 2 


= % f [K o (r/n) + £ K 1 ( r A)]l 0 ( t /^)tp(t)dt 
4 *sO 

R 1 

+ % f K (t/n)[l Q (r/ (j,) -t I 1 (r/ i i)]tp(t)dt 

4 */r 

a 2 

+ ^| ^/ K Q (t/p)[l o (r/n) -£ ^(r/iijjtp^jdt -C 1 p 1 (r) 


(4.8a) 


and 


* 2(r,|J, ) =% / [KlrWtkrt)]! (t/|x)tp(t)dt 
5r 2 V? J Q L ° r 1 Jo 

+ /[^AO + £ K 1 (r/p)]l o (t/n)tp(t)dt 


+ % f K (t/V)[l (r/n) -£ I 1 (r/ Ml )]tp(t)dt -C 2 p 2 (r) 

M* ./r 

(4.8b) 

Now substituting from Eqs. (4.6a), (4.7a), and (4.8a), we form the 
equation 

+ i atl(r,U - ) - 4 

6 r 2 r £>r 4 ■*- 

r R x 

= --§■ y [K Q (r/ w ) K 1 (r/p)]l o (t/fi)tp(t)dt K o (t/n)[l Q (r/n) 

R 2 

- £ I^r/i^MtJdt+Hf / K 0 (t/n)[l 0 (r/n) -£ I 1 (r/n)]tp(t)dt 
r 

-C 1 p 1 (r)-^J- j K 1 (r/^)l o (t/ M ,)tp(t)dt j K o (t/n)l 1 (r/ M ,)tp(t)dt 

& 


r 4 


/ K n (t/n)l (rAOtp(t)dt -% /* K (r/(j.)l (t/u)tp(t)dt - 
rL ° 4V 0 ° ° 
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C 

T 


1 

f K (t/ |i )l o (r/|j 1 )tp(t)dt 

*/r 

; 2 

y K (t/n)l o (r/ l4 )tp(t)dt. 


,,s ./ "o 

R 


(4.9a) 


Upon cancellation of like terms in the right hand side of Eq. (4.9a) 
and the substitution of Eq- (4.5a) for p^(r), Eq. (4.9a) becomes 


6 2 §i(r,|j.) 

6r* 


+ 


r 5r ^ 



M r V) 

n 




(4.9b) 


which is identical in form to Eq. (3-10d). Now substituting from 
Eqs. (4.6b), (4.7b), and (4.8b), we form the equation 


& 2 $ 2 (r,M-) i &§ 2 (r,^) i 


br 


2 + r 


& r u 


R r 

r/ji) +7 K 1 (r/n)Jl o (t/^)tp(t)dt /[V^> 

/ h 2 

K 0 (t/^)[l o (r/n) 

R 1 

“ r I i( r /M')] t P( t ) dt “ c 2P E ( r ) - J K 1 (r/|i)l o (t/(i)tp(t)dt 


r M> 


r R 2 
J K^r/pOl (t/n)tp(t)dt+|j f K (t/n)l 1 (r/n)tp(t)dt 
■R_ ^ Jr 


-% f K 0 (r/n)l o (t/n)tp(t)dt f K (r/n)l (t/n)tp(t)dt 
^ ^ R, 


/*2 

'3y K 0 (t/h) X o( r /4)tp(t)dt. 


(4.10a) 
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Upon cancellation of like terras in the right hand side of Eq. (4.10a) 
and the substitution of Eq. (4.5b) for p 2 (r), Eq. (4.10a) becomes 



(4.10b) 


which is identical in form to Eqs. (4.9b) and (3»10d). 

Having established that $ (r^) and $ 2 (r,p/) obey the same integro- 
differential equation that arises for $(r,p,) in the bare core solution, 
we now wish to determine the boundary conditions satisfied by §^(r,n.) 
and $ 2 (r,g,). Inspection of $^(r,fi) as defined in Eq. (4.6a) at the 
core centerline, r=0, with the consideration that I Q (0) =1, and that 
the limit of tK (t/p,) =0 as t -♦ 0, establishes the boundary condition: 

* (0,u) finite. (4.11) 

Next we evaluate ^(r,^) and $ 2 (r,jj) as defined in Eqs. (4.6a,b) 
at the core-reflector interface, r=R^, obtaining 


^(R^H) = C 1 £ K 0 (R 1 /M.)l 0 (t/n)tp(t)dt +C 2 J K o (t/^)l o (R 1 / M .)tp(t)dt, 


(4.12a) 


§ 2 (R 1 ,p,) = C ± / K 0 (R 1 /n)l 0 (t/n)tp(t)dt +C 2 f K o (t/ M ,)l o (R 1 /| i )tp(t)dt. 

J ° \ 

(4.12b) 

Comparison of Eqs. (4.12 a and b) yields one interface boundary condition as 


^l(®l>M') — 


(^•13) 
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A second interface boundary condition is established by evaluating 
the first derivative terms, Eqs. (4.7a and b) at r=R^, obtaining 


b$i(r,|J.) _ _ c 

r = R 


br 


R 1 % 

/K 1 (R 1 /,)l 0 (t/,)tp(t)dt + ^ J K 0 (t/4)l 1 (R 1 /fi)tp(t)dt 


and 

b#2(r,|i) 

dr 


r-R, 


h c / 2 

J K 1 (R 1 /n)l 0 (t/^)tp(t)dt+^ £ K o (t/n)l 1 (R 1 . 


(4.14a) 

,/n)tp(t)dt 
(4.14b) 


Comparison of Eqs. (4.l4a & b) yields the second interface boundary 
condition as 


br 


b$2 (r,u) 


br 


r=R, 


r =f 


(4-15) 


It is interesting to note that the result of operating on boundary 
condition (4.13) by Eqs. (4.5a, b), thereby converting to the neutron 
density, is the familiar continuity of scalar flux across the Interface 
boundary which is used in neutron diffusion theory solutions. The 
result of performing the same operation on boundary condition (4.15) 
is equivalent to the continuity of the neutron current. This equiva- 
lence is dependent on identical diffusion coefficients in each region 
which holds here due to the identical total cross sections and iso- 
tropic scattering in the postulated two-region model. 

A fourth boundary condition is sought at the outer reflector 
boundary, r=R 2> First we evaluate $ 2 (r,p,), Eq. (4.6b), and its first 
derivative, Eq. (4.7b), at the outer boundary, obtaining 
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$ 2 (R2,uO = C 


and 

b$ 2 ( r ' n) 
5r 




K o (R 2 /fi)l o (t/n)tp(t)dt, 

(4.l6a) 


2a 

4 


J K 1 (R 2 /|i)l 0 (t/n)tp(t)dt. 

1 (4.l6b) 


Inspection of Eqs. (4.l6a,b) as combined below establishes the outer 
boundary condition as 


k(r 2 A*) 

o 2 5r 


+ *( 1 * 2 , 4 ) 


= 0 , 


r = R, 


(4.17) 


which corresponds to Eq. (3>lld) in the bare core solution. 

In summary, we have developed two pseudo neutron distribution 
functions given by Eqs- (4.6a,b) and have shown that they obey the 
same integro-differential equation, Eqs. (4.9b), (4.10b). Also we 
have established four boundary conditions, Eqs. (4.11), (4.13), (4.15), 
(4.17), satisfied by the pseudo neutron distribution functions. 


4.3 Singular Eigenfunction Expansion and Solution 

Since it was shown in the previous section that both $^(r,p,) 
and $ 2 (r,p,) obey integro-differential equations of the same form as 
Eq. (3-10d) in the bare core solution, the general solution for 
$^(r,p.) and $ 2 (r,p) corresponds to the form of Eq. (3*21) which we 
write here as 

v r ’“> - ['Vof-v'v + 

1 

+ P ,v K 0 (r/v) ]“\y (,i)d '' I 4 ' 18 ' 


I =1,2 
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where T| o ^(n,) and T)^(|i) are given by Eqs. (3*19) afld (3*20), respectively. 
Two additional aspects of the bare core solution are present here for the 
pseudo neutron distribution function in the core region. One is that v ^ 
will again be imaginary, giving rise to regular Bessel functions, and the 
second is that the centerline boundary condition, Eq. (4.11), is identical 
to the bare core centerline boundary condition,. Eq. (3.11a), Therefore, 
Eq. (4.18) for ^(r,^) reduces to the form of Eq. (3.22) which we write 
here as 


= b 0 J 0 ( r /| v ol!^ 2T| ol^^ + jf B Ml Q ( r AO^\ v GO d v . (4.19) 


Adopting more specific symbols for the coefficients in the reflector 
region, we write Eq. (4.18) for ( r ^ M- ) as 


» 2 (r,„) . [a o I c (r/v o2 ) + V'o (r/ V>F V>*> 

1 

^ 2 T1 2v ( | x)dv. (4.20) 

The coefficients in Eqs. (4.19) and (4.20) will be obtained from 
the following iterative procedure: 

1. B(v) is initially taken to be zero. 

2. b^ is an arbitrary constant corresponding 
indirectly to the power level which we set 
equal to unity. 

At the core-reflector interface we obtain 


J [A(v)l 0 (r/v) +D(v)K o (r/v)] 


3* a Q as a function of discrete terms and B(v)« 

4. d as a function of discrete terms, a , and B(v). 
o o ' ' 

5* D(v) as a function of discrete terms and B(v)» 
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At the outer boundary we obtain 

6. an inner iterative sequence for A(v) depending 
on a Q , d Q , and D(v)- This sequence provides 
the criticality condition. 

At the core-reflector interface we obtain 

7 . B(v) as a function of a^j d^, A(v)> and D(v)j 
completing the iterative sequence. 

To effect steps 3 through 7 of the iterative solution, extensive 
use will be made of the orthogonality and normalization relations for 
the pseudo eigenfunctions developed in Appendix A. Another useful 
equation which relates the continuum eigenfunctions in different 
regions is written as 


C. 


V* 5 = cT + 




(4.21a) 


To verify Eq. (4.21a), we substitute Eq. (3-20) for 7)^ v ( M- ) obtaining 


T). (n) = -yV 8 . f ; + 1 -C.P f 

Ov ' c i ) 1 v B -n“ 1 J Q 




s(v - n); 


c. 

+ 8(v - n) -jr- &( v - n) 
i 


(4.21b) 


where we have substituted Eq- (2.15) for ^(v). Canceling the like 
terms on the right hand side of Eq. (4.21b) we have 


T|. <(X) = C.P — 3— — -g- 


+ hj (v)S(v - p), 


(4.21c) 


which establishes Eq. (4.21a). 
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Proceeding with the solution, we substitute Eqs. (4.19) and (4.20) 
into the interface boundary condition Eq. (4.15), obtaining 


Vo ( VivnA>) + f B(v)i 0 (V v )^ 2 % v ^) dv 
■ [Vo'Vv’ ,i o‘« ( Vv ) ]* ! V ( ' ) 


+ 


1 

f\ [a(v)I o (R 1 /v) +D(v)K o (R 1 /v)] u 2 T) 2v (^) | iv. 


(4.22) 


Next we operate on Eq. (4.22) with the integral of T^^) over p, 
to yield 


1 

Vo ( V!vl ) f \i^V^ 

'•'o 

1 1 

- /“MW *) £% + ( 5a cr 1 ) 

1 

= [ a o I o (R l /v o2 ) +d o K o ( VV ) ]jA 2 ' ri o2^ ) ^o2 ( ^^ 

1 1 
+ /dv[A(v)l 0 (R 1 /v) +D(v)K q (R 1 /v)] /l\ (p)Tl 2 (^)dfi, (4.25) 

•'u +'Q 


where we have substituted Eq. (4.21a) for T|^ (|j.) - Now we apply the 
orthogonality and normalization relations, Eqs. (A .6), (A. 7) and 
(A* 12), to reduce Eq. (4.23) to the form: 


Vo< E l / M> N ol2 K*(v)ljV 

= [ a oIo( R l/' , o2) + V 0 (B 1 /v o2 )]K o2i 
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At this point we accomplish step 4 of the iterative procedure by solving 
Eq. (4.24) for d Q in terms of the other variables, resulting in 


d ° = KotR./vo [V o'V I ''ol I )"ol2 ‘ Vo< VV> M o2 



(4.25) 


To obtain another equation so as to isolate a , we turn to the second 
interface boundary condition Eq. (4.15) and substitute Eqs. (4.19) 
and (4.20), obtaining 



VVKii^oi^ 



W v)^i 2 Tl lv (fi)dv 


h <V < v 02) - ^ % ( V v )] ) 

1 

+ /[^ I^/v) K^Rj/v)] ,\>)dv. 


(4.26) 


Next we operate on Eq. (4.26) with the integral of over (j, to 

yield 

7^7] J i( R i/l v oil 

1 1 

+ dv ^2\>^ + ( 2 C ; ^ 1 l 0 2^ M ‘) d ^ 

1 

= fc Wv* -^7 
1 1 

+ /dvf^M I^/v) (4.27) 

•Sq L 
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where we have substituted Eq- (4.21a) for T|^ v (n,). Now we apply the 
orthogonality and normalization relations, Eqs. (A. 6), (A.7)> and (A-12), 
to reduce Eq- (4.27) to the form: 




v\ 2 (v) 



(4.28) 


Now, with Eqs. (4.24) and (4.28), we algebraically eliminate d Q and, 
after substituting the Wronskian for the modified Bessel functions, 
obtain the equation for a^ required in step 3 as 



(4.29) 


To develop the equation for D(v) as required in step we return 
to Eq. (4.22), operate with the integral of T^»(|j,) over \i, and thereby 
obtain 



(li)diJ, (4.30) 
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where on the left hand side we have substituted Eq. (4.21a) for T^iG*). 
Now we apply the orthogonality and normalization relations, Eqs. (A. 6) and 
(A. 21), to reduce Eq. (4.30) to the form: 


Vo ( Vlv 

1 

+ /^b( v )i o (r 1 /v)(H^.)v' 2 H 1 v (v') 

= [ACv’J^d^/v') + d(v , )k o (r 1 /v 1 ^n 2 (v’)- 



Now returning to Eq. (4.26), we again operate with the integral of 
, (p,) over |x to obtain 


+ (^^■]s(v' -n)W 


+ -n»dp. 


1 1 

+ A ^ r i ( V v ) /A>>j£ 

o of 1 

= h'Vv’ 

/* 1 

+ f i i^v v) - 


(4.32) 


where on the left hand side we have substituted Eq. (4- 21a) for 

^2v f ' N ° W We apply orthogonality an<i normalization relations, 

Eqs. (A. 6) and (A. 21) to reduce Eq. (4.32) to the form: 
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W'-') +2 ^ 

■\ v (v') 

(R 1 /v , )]»2(v’)- 


I i<V v ') 5f k i ( v,) 


(4.33) 


Now, with Eqs. (4.31) and. (4.33) we algebraically eliminate A(v) and, 
after substituting the Wronskian for the modified Bessel functions, 
obtain the equation for D(v) required in step 5 as 


D(v’) 




1)1, (H,/v t ) 


+ £i 




r i n (R,/v)I 1 (R,/v') 


. I, (R,/*)UK,/v')] j 


(U.3U) 


Substituting Eq. (3*20) for Tj^(v’) in Eq. (4.34) and evaluating the 
delta function leaves the indeterminate form 


F(v,v r ) 


... ^ [” 1q (Ri/v)lj (R 1 / v’ ) I, (R, / v)l n (R-i /v 1 ) 

V 2 - v ' 2 L V 


The function is evaluated by L'Hospital's rule to be 


F(v,v) [l 0 (R 1 /v)l o (R 1 /v) -I 1 (R 1 /v)r i (R 1 /v)]. 


(4.36) 
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To accomplish step 6 of the iterative procedure, we turn to the 
outer boundary condition, Eq. (4.17), and substitute Eq. (4.20) along 
with its derivative for the pseudo neutron distribution function in 
the reflector to obtain 



a a T 0 
o r 'o2 


oo 





iD(v) 4 2 11 2v (n)dv = 0. 


(4.36) 


Separation of Eq. (4*36) into singular and nonsingular parts follows 
directly as in Eqs. (3*23) through (3.27c) in the bare core solution. 

The result of this procedure on Eq. (4.36) is 

1 

J A’(v)Tfe (n)dv = $>), (4.37) 

o 

where 
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A ’ (v) = A(v)l (r/v), and H(v,|i) is defined in Eq. (3-25b). To demon- 
strate that $'(p,) is indeed nonsingular, we further reduce the integral 
containing D(v) by substituting Eq- (3-20) for T^^(|u,) and immediately 
evaluate the delta function to obtain 




2 , 

v dv- 


Next we extract the indeterminate portion of the above integral and 
define it as 


F(v,n) 


1 f K 1 (R a /^)K n (R Q /v) _ K a (R 2 Ai)Kj (Rg/u.) ' 
V s - |X S L V V J 


(4-38) 


Finally, we apply L’Hospital T s rule to Eq- (4-38) for F(v,p,) to obtain 


F(v,v) = Ijr [^(Bg/vJ^CRg/v) - k o (r 2 /v)k q (r 2 /v)] - (4.39) 

With the above demonstration and recalling that H(v,v), Eq. (3*25c), 
is a bounded finite function, we have successfully demonstrated that 
$’(^), as defined above, is a nonsingular function. 

Now the inner iterative sequence for A' (v) required in step 6 
can be constructed in a manner closely analogous to that of the bare 
core solution, Eqs. (3.27a) through (3.29). The associated discussion 
concerning the completeness theorem proof giving rise to a free expres- 
sion for A' (v) and the criticality condition carries over directly. 
Successive approximations for A’(v) are constructed by the following 
iteration scheme : 

A’(v) = 0 
o 


(4.40a) 
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(4.40c) 


(4.40d) 


(4.4oe) 


Eqs. (4.40a) through (4.40d) represent the inner iterative sequence 
for A'(v). The criticality condition, Eq. (4.40e), is checked after 
convergence of the outer iterative sequence for a Q , d , D(v) and B(v) 
and after a simultaneous convergence of the inner iterative sequence for 
A(v)» Among the solutions for two-region problems, this solution is 
unusual in that the criticality condition is developed at the outer 
reflector boundary. 

The remaining portion of the solution is to develop the equation 
for B(v) as required in step 7 . Toward this end, we return to Eq. (4.22) 
and operate with the integral of t (jj. ) over p, to obtain 
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i Y r 2 . 

Vo (R l / l V oll ) /A>1 ( ^ )T W ( ^ + J dvB(v)l 0 (R L /v) U 

- [‘o I o ( VV )td o K o ( V , o2 ) ]/ |,2 V W {| (|j) + ( £a ^ 21 -) 6( ' , '- K) } dw 


-ydv[A(v)l 0 (R 1 /v) +D(v)K o (R 1 /v)]y^ 2 T^ v (n) ^.(p-) 


(c, „ 


S 

(4.41) 


where we have substituted Eq. (4.21a) for , (p,) on the right hand 
side of Eq. (4.4l). Now we apply the orthogonality and normalization 
relations, Eqs. (A .6) and (A. 21) to reduce Eq. (4-31) to the form: 

Bcviyvvjyv) = [VofVv^Vc'V^^Sr 1 -) 

+ [a^'UJRj/v') +d(v , )k o (r 1 /v')]^- n 2 (v') 

+ Ji Lv[a(v)I 0 (R 1 /v) +D(v)K o (R 1 /v)] v’\ v (V). (4.42) 

Rearrangement of Eq. (4.42) yields the required equation for B(v) as 

= '^(r^'vMn^v*) I [Vo ( Vv )+d o K o (R i / V ) ]( ££ cT i ) v, 2 V v,) 

+ [a(v’)I o (R 1 /v’) +D( v')K o (R 1 /v’)]^- N 2 (v’) 

+ ^v[a(v)I o (R 1 /v)+D(v)K o (R 1 /v)]^^ 1^\ V’\ v (v’)i. (4.43) 


We note that the integral term in Eq. (4.43) leads to a singularity. 
The numerical evaluation of such singular integrals is discussed in 


Chapter V. 
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Finally, we identify the various equations which appear in the 
iterative solution for the two-region problem previously outlined. 


1 . 

B(v) is initially taken to be zero. 


2. 

b^ is set equal to unity. 


3. 

a Q from Eq. (4.29). 


4. 

d Q from Eq. (4.25). 


5. 

D(v) from Eq. (4.34). 


6. 

Inner iterative sequence for A(v), Eqs. 

(4.40a) 


through (4.40d), criticality condition 

Eq. (4.40e) . 

7- 

B(v) from Eq. (4.43). 



Upon finding the critical core radius for given values of C^, 

C , and the reflector thickness, the expansion coefficients can be used 
to calculate the neutron density distribution in each region from 
the equations: 


x 

p l( r ) = b 0 J 0 ( r / l v oll^ + >( /B(v)l 0 (r/v)dv 


1 

p 2 ( r ) = Vo< r /V )+d o K o (r/ V } + /[ A ’( v ) + D Cv) K o C r / v)] dv, 


R^riR 2 


(4.46) 


which we have obtained by operating on Eqs. (4.19) and (4.20) as indicated 
in Eqs. (4.5a and b) and applying the normalization of the pseudo 
eigenfunctions as in Eq. (3-l6). 

This completes the solution for the radially reflected critical 
cylinder. Numerical results are presented in Chapter V. 
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CHAPTER V 

NUMERICAL METHODS AND RESULTS 
5«1 Numerical Methods 

In obtaining numerical results from the iterative solutions developed 
in Chapters III and IV, various parameters, functions, and operations 
appearing in the solutions were evaluated by the following methods: 

1. The dispersion function for the discrete eigenvalues, which is 
the transcendental equation, Eq. (2.13), was solved by Newton's 
method . 

2. The regular Bessel functions were calculated with a recurrence 
relationship. 

3. The modified Bessel functions were calculated with a series 
expansion for the smaller arguments and by an asymptotic series 
for the larger arguments. 

4. Singular integrals were treated by subtracting out the 
singularity, evaluating the resulting integral term by Gauss- 
Legendre quadrature and evaluating the derivative term by 
Lagrange interpolation. 

5. Nonsingular integrals were evaluated by Gauss-Legendre 
quadrature . 

Although these methods are generally well-known, for completeness 
we shall give a brief description of each method and make appropriate 
references . 

For c > 1 the discrete eigenvalue as given by Eq. (2.13) is purely 
imaginary. The solutions developed in Chapters III and IV are formulated 
in terms of the magnitude of this imaginary eigenvalue. Therefore, we 
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use Newton ' s iteration scheme as 


v l i+ l= i 


MM ± ) 


(i— 0 ^ 1 , 2 , . . .) 

where A(|v|) = l-c|v| tan "*'(l/|v|) 

and A ' (| vj ) = c|jv| / ( | v ] 2 +l) - tan - 1 (l/|v ] )j . 


(5.1) 


For c < 1 the discrete eigenvalue is real and we apply 
A(v ± ) 

v i+l = v i A «(v i ) > 


Newton's method as 
(5.2) 


(i« 0 ,l, 2 ,...) 

where A (v) = 1 - in^~^ 

and A'(v) = | |~2v/(v 2 -l) - in(-~“)J . 

Initial estimates for | vb | are taken from Fig. 4. The iterative sequences, 
Eqs. (5>l) and (5*2), were run until the error in | \b | was less than 10~^. 

The recurrence relation used to calculate the regular Bessel 
functions is taken from Goldstein and Thaler. 1 ^ The recurrence relation is 


F n + l (x) + p n-l« = ^ V x >- 
The desired Bessel function is 


M-2 

where a = Fo (x ) + 2 X F_ (x ) 

2m 


(50) 

(5.4) 

(5.5) 


and M is initialized at Mo . 

Mo is the greater of M^ and Mg where 

M = if x < 5 

|[l.4x + 6 o/x] if x s 5 
Mg = [n + x/4 + 2 ] 


( 5 . 6 a) 


and 


(5*6b) 
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Figure 4. - Variation of the discrete eigenvalue for 0 <c < 2. (Mean number of second- 
aries per collision). 
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Sequential values F^ '* * * ' ^ 2 * *1* ^0 0X6 eva ^- ua ^ e ^ us i n g Eq. 

(5*3) with F„ = 0 and F„, , = 10~ 5 °. Values of a and J (x) are then 
M M-l n 

computed using Eqs. ( 5 . 5 ) and (5.4), respectively. The computation is 

repeated for M + 3 and the values of J^(x) for M and M + 3 are compared. 

If these values differ by less than 10 the value of J n ( x ) i- s accepted. 

If the test fails, the computation is continued by adding 3 to M and 

using this as a new value of M. 

The methods for calculating the modified Bessel functions are 

12 

taken from Abramowitz and Stegun. For values of the argument less 
than 8 . 5 , we use the ascending series: 


1 ( z ) = 1 + £ L 

° wf 


(^ 2 )\ ( i * 2 ) 5 




(3!)- 


( 5 . 7 ) 


and 

1 z 2 

K (z) = -Un(£z) + v3 I ( z ) + JL 

0 0 TITT 

+ ^ 1+ 2^ z2 ) 2 + ( 1+ |+ t ) (7* z2 ) 5 + . . (5*8) 

T217 2 TTTy 3 

where y is Euler 's constant* The first order modified Bessel functions 
are calculated from the differentiated forms of Eqs* (5*7) and (5*8) 
according to 

i-l ( z ) = i; ( z ), ( 5 . 9 a) 


and (z) = -K^ (z). 


(5-9b) 


The expansions were carried out to thirty terms or until a major under 
flow results from the calculation. For values of the argument greater 
than 8 . 5 , we use the asymptotic series 
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(p-l)(p-9)(p-25) 

31 (8z)3 + ' ’ * 

and 


v>*V£ + 


p 

where to. = 4n . 


+ (P-D(P-9)(p-25) + 

37 (8z)3 


The series were carried out to thirty-two terms. 


(5.10) 


(5.11) 


Principal value integrals were treated by subtracting out the 

14 


singularity as shown in Metcalf and Zweifel/ 


'/w -i 


[f(h') - F(n)i a* 

p’ - P 


+ F(p) £n 


13L 

P 


and approximating the principal value by 
1 M 

1 ~ y w 


-/ 


p '-P 




F(p .)-F(p.)' 


^ d"**i 


W ± dF(M') 


dp 


P=Pi 


(5.12) 


F(u 1 )jgn 


. P-i 


(5.13) 


where the singularity is at p^ and the Vs and p's are the weight 

factors and abscissas of Gauss -Legendre quadrature adjusted to the 

half-range as shown below. The derivative term in Eq. (5.13) is 

evaluated by Lagrange's differentiation formula: 

M 

F'(^) = E 4 + (5.u) 

k=o 

« M 

where £ * (p.) = £ TT (Pi-P n ) 

j=o n=o (5.15a) 

j^k M ’ 

(P ± -P k ) (Pi-Pj) IT (P k -P n ) 
n-o 
n/k 
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M 


and (n ± ) = TT (m^-M^) 


n— o k , ri M+l 

n ^ ± (M+l ) l djj- 


(5.15b) 


|i=?, 5 ? * H M 

Using M Gauss -Legendre quadrature points, the error in calculating 
F*(p.) by Eq. (5.14) is proportional to the (MKLl th order derivative of 
the function as shown by Eq. (5.15b). At the beginning of the computer 
program, the coefficients t' (|i.) are calculated by Eq. (5. 15a) and stored. 

K X 

This procedure significantly simplifies the subsequent numerical evalu- 
ation of the principal value integrals in the iterative solutions. 

Nonsingular integrals were evaluated by Gauss-Legendre quadrature 
over the half- range as 


/ I M 

F (p,) dp. = S w F(|i, ), 
o 1=1 1 1 


(5.16) 


where = (X^ + l)/2, 

and w . = W./2, 
l r } 

and being the full-range Gauss-Legendre abscissas and weight 

12 

factors listed in Abramowitz and Stegun. 


5.2 Numerical Precision 

The objective in seeking numerical results from the iterative 
solutions of Chapters III and IV was to obtain critical dimensions 
accurate to a minimum of six significant figures. Results of this high 
precision can be used as input data for evaluating design analysis 
programs as is demonstrated in Chapter VI. Therefore, the computer 
programs used to obtain these results were written entirely in double 
precision. The calculations were performed on IBM- 7094-11 and CDC-6600 
computers. A minimum of twelve significant figures were carried in 
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performing the calculations . 

Where possible , comparisons were made between the results of methods 

described in section 5*1 and tabulated values in the literature. The 

discrete eigenvalues agreed with the seven- significant- figure values of 

Kowalska. 1 ^ Sample values of the Bessel functions agreed with the 

12 

tabulated values in Abramowitz and Stegun to at least nine significant 
figures . 

The bare and reflected core calculations were done over Gauss- 
Legendre quadratures of twenty- four and forty points respectively. The 
variation of the critical core radius with the number of Gauss-Legendre 
points is shown in Table I. It indicates that twenty- four Gauss-Legendre 
points are sufficient for seven significant accuracy in the bare core 
radii. The use of forty Gauss-Legendre points provides seven significant 
figure accuracy for the large reflected cores and six significant figure 
accuracy for the smaller configurations. This reflects the increasing 
importance of the continuum contribution and its integrated quantities 
in the solution for the smaller configurations . 

At each trial value in the critical radii searches, the iterated 
quantities in the solutions were converged to nine significant figures 
before the criticality conditions were evaluated. The criticality 

conditions, Eqs. (3*29) and (4.40e), were considered satisfied when 

-9 

the value of the integral was less than 10 . 

3.3 Bare Core Results 

Critical radii calculated with the one-region cylinder solution of 
Chapter III, along with results of other analytic solutions found in the 
literature, are listed in Table II. Using the results of the present 



Case Description 


C 1 

C 2 

R 2 - R 1 

8 

1.02 

Bare 

— 

9.0432542 

2.0 

Bare 

— 

.66862087 

1.1 

0.9 

10.0 

2.795180 

1.4 

.85 

1.0 

1.141910 


of Gauss-Legendre Points 


24 

4 o 

96 

9.0432547 

9.0432547 


.66861281 

. 66861285 


2.795122 

2.795120 

2.795120 

1.141755 

1 

1.141750 

1 . 141748 










TABLE II. - CRITICAL RADII IN MEAN FREE PATHS FOR BARE CYLINDERS 


! c 

i 

! 

Present 

Solution 

Carlson- 

Bell 16 

1.02 

9.043255 

9-0433 

1.05 

5 . 411288 

5 • 4ll8 

1.1 

3.577391 

3.5783 

1.2 

2.287209 

2.2884 

1.4 

1.396979 

1-3973 

1.6 

1.020839 

1.0209 

1.8 

0.807427 

O.8067 

2.0 

O.668613 

0.6673 


Hendry 

F 3 G 8 


5 .4l4 


Herabd 

IT4 


9.04458 


1.39699 

1.02085 

0.80743 

0.66862 
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work as reference values , we see that the values of Carlson and Bell are 

most accurate for the largest system (c — 1*02) and least accurate for 
the smallest system (c = 2.0)# Carlson and Bell used the extrapolated 
endpoint method for the larger systems, R >1-5- For the smaller systems, 
they interpolated between values calculated with the extrapolated endpoint 
and variational methods* The values of Hendry were obtained by Fourier 
expansion of the neutron distribution function in one of the angular 
variables and solution of the resulting equations using integration by 

lQ 

Gauss quadrature. The IT^ method of Hembd is based upon a Fourier 

transformation of the integral equation. It appears to be particularly 

effective, especially for the smaller systems. The reduced accuracy of 

18 

the IT method for the larger systems is recognized by the author. 

The neutron density or total flux was calculated from Eq. (3*3l). 

Table III lists the neutron density relative to the centerline value. 

It is seen that the radial drop-off in the neutron density decreases 

with a reduction in the size of the system. Figure 5 shows the neutron 

density distribution for the smallest case, c = 2.0, calculated by the 

method presented here and by various order discrete ordinate calculations. 

19 

The discrete ordinate calculations were performed with the TDSN program 
using a moment modified quadrature. The neutron density distribution 
from the calculation agrees very closely with the exact distribution. 
However, the excess system multiplication calculated by discrete ordinates 
indicates the increasing importance of the error in the neutron density 
distribution as the discrete quadrature order is reduced. 

The first terra of the neutron density, Eq. (3*31), is the asymptotic 
density arising from the discrete mode. The asymptotic density 



1.02 


1.05 


1.1 


1.2 


r /Bc 

1.02 

0.0 

1.0000 

0-25 

0.9236 

0.50 

0.7118 

0.75 

0.4124 

0.85 

0.2821 

0.91 

0.2033 

0.95 

0.1502 

0.98 

0.1086 

1.0 

0.0747 


1.0000 

1.0000 

0.9299 

0.9381 

0.7340 

0.7561 

0.4522 

0 . 4922 

0.3267 

0.3718 

0.2492 

0.2960 

0.1958 

0.2430 

0.1531 

: 0.1999 

0.1179 

! 0.1641 


1.0000 

0.9433 

0.7819 

0.5399 

0.4263 

0.3535 

0-3016 

0.2589 

0.2233 


1.4 

1.6 

1.8 

1.0000 

1.0000 

1.0000 

0.9508 

0.9550 

0.9578 

0.8093 

0.8248 

0.8352 

0.5918 

0.6218 

0 . 6421 

0.4668 

0.5223 

0.5466 

0 . 4i8i 

0 . 45 66 

0.4830 

0 . 3686 

0 . 4o88 

0.4366 

0.3273 

0.3688 

0.3976 

0.2926 ! 

0.3351 

0.3646 
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/Exact (k - 1. 0000) 
ls 16 (k = 1.0029) 



0 .2 .4 .6 .8 1.0 

Fraction of radius, R/R 0 

Figure 5. - Neutron density 
distribution (C = 2. 0). 
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corresponds to the diffusion theory solution based on the exact diffusion 
coefficient. Table IV lists the ratio of the asymptotic to total 
neutron density as a function of position. The value of this ratio 

at the outer boundary for the case c = 1.02 agrees with the corres- 

3 

ponding value for the Milne problem with c = 1.0. As anticipated, 
the error in the asymptotic density is most severe on the boundary of 
the smallest system. 

.SJt Reflected Core Results 

Critical core radii as a function of core and reflector multiplying 
properties and reflector thickness are presented in Table V. We note 
the anticipated increase in critical core radius with a decrease in 
either reflector multiplication or reflector thickness. 

The critical dimensions of two widely varying cases in Table V 
were used in one-dimensional discrete-ordinates and diffusion theory 
analyses. The results of these analyses are given in Table VI. The 

19 

discrete-ordinates calculations were performed with the TDSN program 
using a moment modified quadrature. For both types of calculations 
it is seen that the variation of the effective multiplication factor 
from critical is strongly dependent on the amount of neutron absorption 
in the core. Thus it appears that in evaluating numerical methods by 
comparison with exact analyses, one should include a realistic amount 
of absorption in choosing the cross sections. A relatively high order 
of angular quadrature is required for good discrete-ordinates analysis 
of the small system. Diffusion theory is adequate for analyzing the 
large system which more closely corresponds to power reactors. 



TABUS IV. - ASYMPTOTIC RELATIVE TO TOTAL NEUTRON DENSITY 







TABLE V -CRITICAL CORE RADII IN MEAN FREE PATHS FOR REFLECTED CYLINDERS 


Case 


Reflector Thickness 

(mfp) 


c i 

c 2 

1 

3 

6 

10 

20 

1.02 

0.99 

8.160095 

6 . 981947 

6. 220461 

5.914592 

5.814384 

1.02 

• 95 

8.286641 

7.621081 

7.431634 

7 . 410246 

7.409152 

1.02 

• 90 

8.411027 

8.036229 

7 . 981761 

7.979325 

7.979288 

1.02 

.85 

8.508960 

8.276755 

8.255960 

8.255474 

8.255471 

1.05 

0.99 

4.618945 

3.772231 

3.325398 

3.159486 

3.105828 

1.05 

• 95 

4.724733 

4.203900 

4.068247 

4.053H7 

4.052340 

1.05 

• 90 

4.831448 

4.520869 

4.477548 

4.475609 

4.475579 

1.05 

-85 

4.917318 

4. 718400 

4.700988 

4.700581 

4.700578 

1.1 

0.99 

2.899832 

2.323758 

2.057837 

1.961971 

1.930929 

1.1 

.95 

2.982381 

2.602577 

2 . 511010 

2 . 500863 

2.500339 

1.1 

• 90 

2.068188 

2.828581 

2.796556 

2.795120 

2.795098 

1.1 

.85 

3.13909 

2.97985 

2 . 96628 

2.96596 

2.96596 

1.2 

0.99 

1.76189 

1.41589 

1.26942 

1.21684 

1.19964 

1.2 

•95 

1.81837 

1.57685 

1 . 52221 

1.51616 

1.51585 

1.2 

■ 90 

1.87916 

1.71823 

1.69762 

1.69670 

1.69668 

1.2 

• 85 

1.93105 

1.81967 

1.81047 

1.81025 

1.81025 

1.4 

1 

0-99 

1 1.03831 

.853831 

.778234 

.750720 

.74i64o 

1.4 

• 95 

1.07179 

.938396 

.909237 

.905989 

.905820 

1.4 

• 90 

1.10898 

1.01652 

1.00504 

1.00452 

1.00451 

1.4 

.85 

1.14175 

1.07560 

1.07028 

I.07015 

1.07015 




TABLE VI -NUMERICAL ANALYSES OF CRITICAL CONFIGURATIONS 


Case Description 

Cl 

. 1.02 

1.02 

1.4 

1.4 

C 2 

■ 99 

.99 

.85 

• 85 

* — 1 

1 

CM 

20.0 

20.0 

1.0 

1.0 

Z al 

.0 

.25 

.0 

.25 

Effective Multiplication Factor* 

Di screte Ordinate s 

s 4 

1.00695 

1.00051 

1.02280 

1.01391 

s 8 

1.00288 

1.00021 

1.00753 

1.00462 

s l6 

1.00224 

1.00016 

1 . 00217 

1.00134 

s 32 



1.00074 

1.00045 

s 64 



1.00035 

1.00021 

Diffusion Theory 


0.98310 

0.99873 

0.80653 

, 0.87137 


w To be compared with = 1. 
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The neutron density distribution as calculated by the present 
solution and by discrete ordinates is presented in Fig. 6. The distri- 
bution from the calculation agrees with the exact values to at 
least three significant figures at all points. The tendency of the 
distribution as calculated by discrete ordinates to peak away from the 
core centerline is most pronounced for the calculation but persists 
even into the calculation. Here again, the excess multiplication 
calculated by discrete ordinates indicates the increasing importance 
of the error in the neutron density distribution as the discrete 
quadrature order is reduced. 

The values of tabulated neutron density distributions for a number 
of cases, as supplementary to the critical dimensions given in Table V 
when applied as analytic standards, are given in Appendix C. 
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Figure 6. - Neutron density distribution. C^, 1.4; 0% 0.85; reflector 
thickness, 1MFP' 



CHAPTER VI 


APPLICATIONS AS AN ANALYTIC STANDARD 
6.1 The Analytical Model 

For this study, a model was chosen which, in the context of the two- 
region solution developed in Chapter IV, most closely approximates the 
Advanced Power Reactor concept described by Whitmarsh^ and shown in 
Fig. 7. The operational requirements for this space power reactor 
concept are that it provide 2.17 thermal megawatts for 50,000 hours with 
a coolant outlet temperature of 1222°K to a Brayton-cycle power conversion 
system. The design requirements of compact size, long core lifetime, 
and high operating temperature dictate the materials appearing in the 
system and therefore its neutronic characteristics. The concept 
employs highly- enriched uranium nitride fuel, tantalum-based alloy clad 
and structural material, lithium coolant, and molybdenum reflectors. 

This material composition leads to a very hard spectrum, fast reactor 
with a median neutron energy of 0.44 MeV. 

Two major problem areas arise in the neutronic analysis of this 

design. The first problem area concerns the adequacy of the neutron 

cross sections for these relatively little-used materials. Estimates 

of the reactivity biases due to the cross sections have been obtained 

20 

from the analysis of small, fast-spectrum critical assemblies containing 

, 21 
these materials and are reported in Mayo and Lantz. Unfortunately, 

it is not possible to separate the error in the reactivity due to the 

cross sections from that which arises from the approximations inherent 

in the method of analysis. This difficulty brings us to the second 

major problem area, the analytical bias. Since the analyses leading to 



Outlet-end, control- 
drum drive nozzles 



Figure 7. - Advanced paver reactor concept 
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the critical dimensions presented in Chapter V are based on hypothetical 
c values, these dimensions are particularly suitable for use in studying 
various contributions to the analytical bias in the design methods. 

At the high neutron energies present in this reactor, the cross 

sections are much smaller than those associated with thermal reactors. 

Thus the average neutron mean free path is of the same order of magnitude 

as the physical dimensions of the system. The resulting relative 

proximity of any position in the system to the core and reflector 

boundaries precludes the use of neutron diffusion theory in performing 

design analysis. Consequently, the design analysis has been performed 

19 22 

with the discrete- ordinates transport theory programs TDSN, ANISN, 

23 

and DOT-IIW. J 

The azimuthal asymmetry of this drum- controlled concept necessitates 
two-dimensional calculations in the analysis of quantities such as 
control swings and radial power distributions. With these large spatial 
descriptions, computer running time and storage considerations have 
restricted the analytical model to the angular quadrature approxi- 
mation. The first investigation of the present study concerns the effects 
of quadrature order and quadrature type upon the reactivity. The second 
investigation is concerned with a consistent discrepancy between the 
multiplication factors for a single configuration calculated with x-y 
and R - 0 geometry descriptions. Each of these investigations employs 
the analytical model described in Table VII. The neutron cross sections 
are approximate one group values. They have been derived from spatial 
and energy group flux-averaged reaction rates taken from a multigroup 
discrete- ordinates calculation of the Advanced Power Reactor concept. 
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Although the cross sections only approximately represent the actual 

system, in the context of the one group model they accurately produce the 

core and reflector c values. That is, neutron absorption and secondary 

neutron production are conserved. To meet the criterion of identical 

neutron mean free paths in the core and reflector regions, the cross 

sections are adjusted to a normalized total macroscopic cross section 

of one cm - '*" in each region. The core and reflector dimensions in mean 

free paths closely correspond to the one-group average values of the 

actual configuration. These dimensions, along with the c values, 

represent one of the cases in Table V. 

6 ■ 2 Discrete-Ordinates Quadrature Study 

This study is concerned with the effects of quadrature type and 

quadrature order as used in calculating the analytical model in R-0 

23 

geometry with the DOT-IIW discrete- ordinates transport program. 

The interest in R-0 geometry arises from the sixty degree azimuthal 

symmetry present in Advanced Power Reactor configuration. The effects 

of quadrature type are studied in the approximation corresponding 

to the approximation level employed in the design analysis. 

In addition to the even moment quadrature set ordinarily used with 

the DOT-IIW program, seven other quadrature sets are included in the 

study. The various prescriptions under which the quadrature points were 

selected are described in the literature. The even moment, odd moment, 

2k 

level moment, and P^ T^ quadrature sets are from Lathrop and Carlson. 

25 

The projection invariant sets A and B are from Carlson. The DP^ 

26 

quadrature set is from Carlson and Lee. And the moment modified 

27 

quadrature set is from Carlson. 
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The quadrature sets are presented in Tables VIII through XV. The 
points are located on the surface of the unit sphere an octant of which 
is shown in Fig. 8 . The tabulated values include eight positive values 
of the direction cosine Tj followed by eight negative values. Thus the 
values include the entire upper half of the surface of the unit sphere. 
The points are located along fixed values of the direction cosine §. 

It should be noted that these R-0 quadrature sets are orthogonal to the 
more familiar X-Y and R-Z quadrature sets which are based upon fixed 
values of the direction cosine called T] in our convention. 

The R-0 geometry used to describe the analytical model in the 
DOT-IIW calculations is shown in Fig, 9 . The largest radial mesh 
intervals are 0.26 cm (MFP) in the core and 0.3 cm in the reflector. 

Much smaller mesh intervals are grouped near the center and at the 
boundaries. The angular mesh interval is one degree. Perfect reflection 
boundary conditions were applied at the top, center, and bottom bounda- 
ries. A vacuum boundary condition was applied at the outer boundary. 

The cross sections used appear in Table VII. 

The results of this study are given in Table XVI. From the view- 
point of quadrature type, it is reassuring to observe that the even 
moment set which is ordinarily used in DOT-IIW and the moment modified 
set which is ordinarily used in TDSN give the best results. 

From the viewpoint of quadrature order, it is somewhat disturbing 
to observe the 0.6$ Ak/k reactivity bias associated with the results 
given by the standard design analysis approximation. However, this 
information can be used to obtain an estimate of the reactivity bias 
associated with the design analysis of the actual configuration. 



Table VIII 


Even Moment, S^, R-0 Quadrature Set 



*The sequence listed here corresponds to the input sequence 
used in the DOT-IIW code. The negative values of T] in the 
second quadrant are positioned in the same order as the 
first quadrant values shown in Fig. 8. 







Table IX 


Level Moment, S 4 , R-0 Quadrature Set 


i 

w i 

1*111,1+8 

**1,1+8 

•1 

0.0 

0.00000001 

-0. 95006840 

2 

0.08333333 

0.31276157 

-0.89711210 

3 

0.08333333 

0.89736271 

-0. 31204180 

k 

0.08333333 

0.89736271 

0.31204180 

5 

0.08333333 

0.31276157 

0.89711210 

6 

0.0 

0.00000001 

-0.44180300 

7 

0.08333333 

0.31276157 

-0.31204180 

. 8 

0.08333333 

0.31276157 

0.31204180 




Table X 


Moment Modified, R -0 Quadrature Set 


i 

“i 



i 

0.0 

0.00000001 

-0.93743690 

2 

0.08333333 

0.35874166 

-0.86607870 

3 

0.08333333 

0.86607872 

-0.35874166 

4 

0.08333333 

0.86607872 

0.35874166 

5 

0.08333333 

0.35874166 

0.86607170 

6 

0.0 

0.08650678 

-0.49236600 

7 

0.08333333 

0.35874166 

-0.34815530 

8 

0.83333333 

0.35874166 

0.34815530 




Table XI 


Odd Moment, S^, R-6 Quadrature Set 


i 


1 ^1,1+0 


i 

0.0 

0.00000001 

- 0.95522640 

2 

0 . 08 J 33333 

0.29587580 

-0.90824830 

3 

0.08333333 

0.90824826 

-0.29587590 

4 

0.08333333 

0.90824826 

0.29587590 

5 

0.08333333 

0.29587590 

0.90324830 

6 

0.0 

0.00000001 

-0. 418431 59 

7 

0.08333333 

0.29587580 

-0.29587590 

8 

0.08333333 

0.29587580 

0.29587590 




Table XII 


P^T^, S^, R-0 Quadrature Set 


i 

“i 

N 1,1+8 

p i,i+8 

i 

0.0 

0.11676916 

-0.94045250 

2 

0.08l5l8l4 

0.54041858 

-0.86884590 

3 

0. 08151814 

0.86096578 

-0.55988780 

U 

0.08151814 

0.86096578 

0.55988780 

5 

0.08151814 

0.54041858 

0.86884590 

6 

0.0 

0.11548774 

-0.50857410 

7 

0.08696569 

0.54041858 

-0.55947480 

8 

0.08696569 

0. 54041858 

0.55947480 




Table XIII 


DP^, S^, R-6 Quadrature Set 


i 

®1 


%i + 8 

i 

0.0 

0.00000001 

-0.97741590 

2 

0.125 

0.57735032 

-0.78867510 

3 

0.0625 

0.95429742 

-0.21132490 

4 

0.0625 

0.95429742 

0.21132490 

5 

0.125 

0.57735032 

0.78867510 

6 

0.0 

0.00000001 

-0.61481020 

7 

0.0625 

0.57735032 

-0.21132490 

8 

0.0625 

0.57735032 

0.21132490 




Table XIV 


Projection Invariant Set A, R-0 Quadrature Set 


i 


N 1,1.8 

^i,i + 8 

i 

0.0 

0.00000001 

-0.94280904 

2 

0.08333333 

0.33333333 

-0.88191710 

3 

0.08333333 

0.88191710 

-0.33333333 

4 

0.08333333 

0.88191710 

0.33333333 

5 

0.08333333 

0.33333333 

0.88191710 

6 

0.0 

0.00000001 

-0. 47140452 

7 

0.0S333333 

0.33333333 

-0.33333333 

8 

0.08333333 

0.33333333 

0.33333333 




Table XV 


Projection Invariant Set B, R-0 Quadrature Set 


i 

0 ). 

1 

Ni,f8 


i 

0.0 

0.00000001 

-0.71007688 

2 

0.03353333 

0.09175171 

-0. 70412415 

3 

0.08333333 

0. 7041241 5 

-O.09175171 

4 

0.03333333 

0.704i24i5 

0.09175171 

5 

0.08333333 

0.09175171 

0. 70412415 

6 

0.0 

0.00000001 

-O.99578192 

7 

0.08333333 

0.70412415 

-0.70412415 

8 

0.08333333 

0.70412415 

0.704i24i5 
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Figure 8. - Direction mesh for $4 quadrature. 
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\ 

u Core- reflector boundary 


Figure 9. - R-0 description of the analytical model. 
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Table XVI 

APR Model Quadrature Study' 


Quadrature 

Type 

Quadrature 

Order 

k ** 

K eff 

DP 1 

it 

1 . 008 l 

P 3 T 4 

it 

1.0065 

Level Moment 

it 

1.0095 

2 

Even Moment 


1.0062 

Even Moment 

8 

1.0016 

Odd Moment 

it 

1.0112 

3 

Moment Modified 

it 

1.0064 

Pro j . Inv. Set A 

4 

1.0074 

Proj. Inv. Set B 

Ip 

1.0211 


All calculations performed with the DOT-IIW 
program in R-(9 geometry using the diamond 
difference model. 

P 

Normally used with DOT-IIW. 

"3 

J Normally used with TDSN. 


eff 




■r 


* To be compared with k 


= 1 . 
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The difference between the multiplication factors given by the and 
Sq calculations in Table XVI is 0.46$ Ak/k. The corresponding value 
given by the analysis of the actual configuration is 0.5$ Ak/k. The 
larger value given by the analysis of the analytical model reflects 
the greater importance of the radial leakage in the axially infinite 
model. By scaling the 0.6$ Ak/k reactivity bias by the ratio of these 
differences, we obtain an approximate value of 0.4$ Ak/k reactivity 
bias to be used with the standard S^, R-9 design approximation. This 
value is significant when compared to the 0.52$ Ak/k total reactivity 
requirement for the first 10,000 hours of operation as calculated by 
Whitmarsh. 10 

6.5 X-Y Boundary Approximation Study 

The motivation for this study is the observation that, for the 

type of fast reactors under design consideration, an R-9 calculation 

of the multiplication factor consistently produces a higher value than 

does an X-Y calculation of the same parameter. This phenomenon was 

21 

noted by Mayo and Lantz in their pre-critical analyses. In the 
neutronic design analysis of the Advanced Power Reactor concept by 
Whitmarsh, 10 this discrepancy is determined to be worth 0.65$ Ak/k 
in reactivity. 

An X-Y representation of a circular boundary is done in a stepwise 
manner. This results in an outward displacement of the core material 
into the reflector region accompanied by an inward displacement of the 
reflector material into the core region. Also the length of a stepwise 
boundary is a factor of 4/tt greater than the actual boundary. The 
combined effects of having the fuel material in a lower importance 
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region and having a higher probability of neutron leakage from the core 
will result in a negative reactivity effect . This will be somewhat off- 
set by the positive reactivity due to the possible unrealistic first- 
flight neutron reentry into the core. 

On the other hand, an R-0 geometry is capable of describing 
circular boundaries exactly. Therefore it might be assumed that the 
lower multiplication factors calculated in X-Y geometry are in error 
due to the effects discussed above. However, without an absolute 
standard it is not possible to assess the accuracy of either geometrical 
representation. This demonstrates the usefulness of applying the highly 
precise analytic values of Chapter V in resolving this discrepancy. 

The approach taken in performing the study was to calculate the 
analytical model in systematically finer X-Y geometrical representations 
and note the improvement in the multiplication factor due to the 
reduction of the effects discussed above. A computer program was 
written to determine the X-Y boundaries and provide the input for the 
discrete -ordinates program, DOT-IIW. 

The program selects the X-Y boundaries according to three criteria. 

1. The areas of the core and reflector regions are preserved. 

2. Each of the interstitial areas formed by the X-Y and circular 
boundaries is no greater than a maximum input value. 

3. Paired sets of successive interstitial areas exterior and 
interior to the circular boundary are the same size. 

These criteria are effected in the following manner. First, we consider 
two successive interstitial areas as shown on the next page: 
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to the input value of the maximum interstitial area. The value of Xg 


is found by integration of 
x- 


A 

obtaining 


=/ x * (^i - V r2 - x2 ) dx > 


(6.1a) 




X 2 ,~2 2,# , R 2 


sin ' 1 (r) - 


" y l X l + ~ < R " x l^ 2 + J- sin , (6.1b) 

and solving the transcendental Eq. (6.1b) for Xg by Newton's method. 

Having found x.^, w® determine ^2 by 72 = V ^ . 


Next, the value of y^ is found by integration of 


■f 72 ( Vr 2 -7 2 -x 2 )dy, 


(6.2a) 


obtaining 


A = 


t* 2 -#* + f S1B ' 1 (r) - wz 

and solving the transcendental Eq. (6.2b) for y^ by Newton's method. 


(6.2b) 
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This procedure is repeated until the unknown x is found to he 

greater than R/V "2 or the unknown y is found to be less than R /J~2~ . 

The point, x = y = R// 2 , corresponds to where the forty- five degree 

line drawn from the center intersects the circle. The symmetry about 

this line allows us to reflect the boundaries found for the upper 

octant of the circle into the lower octant. That is, x = y. , y = x , 

3 n 1 3 n 1* 

x^_^ ~ ^2 3 ^n 1 = X 2 , e ^ c * -^e ori -*-y P r °t>lem remaining is the approach 

to the interface point. If this point is passed on an x search, x^ 

is determined by 

r 2 ["2 sin" 1 x i - n/2] , (6.3) 

X 2 = 2(x i -y i ) L R - J 

which is found by solving Eqs. (6.1b) and (6.2b) simultaneously and 
substituting 

x i = y 3 an<1 ^ = y 2 . 

A similar procedure cannot be used upon passing the interface 
point on a y search, since criterion 1 could then not be satisfied. 

In this instance, the previous x search calculation is re-performed 
on successively finer values of the interstitial area until the 
subsequent y search produces a value greater than R // 2 . Then 

the final value of is found by Eq. (6.3). 

Two of the four x-y descriptions of the analytical model which 
have been determined by the above procedure are shown in Figures 10 
and 11. The input criterion for the maximum interstitial area was 
0.04 MFP squared for the coarsest mesh. This value was successively 
halved to a value of 0.005 MFP squared for the finest mesh. 
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Figure 10. - X-Y description of the analytical model, 16 x 16 
mesh. 
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1.84 

1.92 


Figure 11. - X-Y description of the analytical model, 22 x 22 
mesh. 
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The calculations were performed with the DOT-IIW discrete-ordinates 
program using the diamond difference model and the even moment, S^, 
quadrature. Thus, these calculations differ from the even moment, S^, 

R-0 calculation of Table XVI only in the geometry description. 

The results of the x-y calculations are given in Table XVII. As 
anticipated, the calculated value of the multiplication factor is low 
and it improves as the x-y mesh is refined. However, it was unanticipated 
that even the least accurate x-y calculation yields a better value of 
the multiplication factor than that given by the R-0 calculation. 
Furthermore, the results of the finer x-y mesh calculations appear 
to be approaching 0-1 io Ak/k below critical. This value is better than 
that given by the Sg, R-0 calculation. The general conclusion of this 
study is that the source of the discrepancy between the results of x -y 
and R-0 calculations lies in the inaccuracy of the R-0 calculation. 
Therefore, it is recommended that design calculations fran which highly 
precise values of the multiplication factor are required, be done in the 
x-y geometry. It is noted that it was not possible to resolve the 
source of this discrepancy without the highly precise results of the 
analytical solution presented in Chapter IV. 

As shown in Chapter V, both the error in the multiplication factor 
and the error in the neutron density distribution from R-0 calculations 
are quadrature dependent. Therefore, it is reasonable to assume that 
these errors are related and that possible methods to eliminate the 
"centerline dip" and improve the neutron density distribution as 
calculated by low quadrature order, discrete-ordinates calculations in 
R-0 geometry would be a worthwhile subject for further research. 



Table XVII 


X-Y 


Step Boundary Approximation Study 


1 


2 

Interstitial Area (MFP) 

Mesh Size 

k * 

eff 

0.04 

16 X 16 S 

0.99539 

0.02 

22 X 22 

0.99718 

0.01 

28 X 28 

O .99852 

0.005 

58 X 58 

0.99895 


1 A11 calculations performed with the DOT-IIW program 
using an even moment, S^ quadrature and the diamond- 
difference model. 
tt To be compared with = 1 . 




CHAPTER VII 


CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 

The conclusions drawn from this study fall into two general 
categories. The first category concerns the successful development of 
the singular eigenfunction solution for the reflected critical cylinder 
through the use of new relations for the singular eigenfunctions. The 
second category concerns the application of the highly precise numerical 
results obtained from this solution in assessing the accuracy of more 
approximate methods of transport analysis. 

In evaluating the singular eigenfunction solution for the bare 
critical cylinder by Mitsis, it was found that his solution contained 
an unbounded function. This function was reformulated and the continuum 
expansion coefficient was redefined to obtain a converging solution. 
Furthermore, it was shown that the procedure taken by Mitsis in decompos- 
ing the singular eigenfunctions into the sum of the Case plane geometry 
eigenfunctions was unneccessary and that it introduced additional 
numerical complexity into the solution. This procedure was avoided 
through the use of a completeness proof for the singular eigenfunctions 
which naturally arise in the solution of the transformed integral 
transport equation in cylindrical geometry. 

In developing the singular eigenfunction solution for the reflected 
critical cylinder, application was made of the reformulated function and 
the completeness proof mentioned above. Also, extensive use was made 
of orthogonality and normalization relations found for the singular 
eigenfunctions. The solution is in the form of two coupled iterative 
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sequences with which the expansion coefficients and ultimately the 
critical dimensions can be found to any desired degree of accuracy. 

Critical core radii accurate to seven significant figures are 
presented for bare cylinders and compared with the results of other 
analytic methods. It was found that for small systems, C >1.6, 
values that had previously served as analytic standards are accurate 
to only three significant figures. 

Highly precise values of critical core radii as a function of 
core and reflector multiplying properties and reflector thickness are 
presented. Comparison of these results for two widely varying systems 
with the results of numerical solutions demonstrates the relative 
accuracy of various design approximations. High order discrete- 
ordinates calculations gave good results for both systems while the 
diffusion theory was inadequate for analyzing the smaller system. It 
was shown that errors in the multiplication factor and the neutron 
density distribution as calculated by discrete ordinates are dependent 
on the order of angular quadrature used. Also, it was shown that the 
introduction of a realistic amount of neutron absorption in making-up 
the core c values significantly improved the results of the numerical 
calculations . 

Additional conclusions were reached in a study of design analysis 
methods using the highly precise results of the present solution as 
analytic standards. These analytic standards are particularly effective 
because of their lack of dependence upon measured values of neutron 
cross sections. Three conclusions concerning the neutronic design 
analysis of the fast reactor concept under consideration were made. 
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First, of eight different types of angular quadrature sets studied, the 
even moment set ordinarily used in the DOT-IIW discrete-ordinates program 
is as effective as any of the other sets studied when applied in the 
approximation in R-6 geometry. Second, using this level of design 
analysis introduces an estimated analytical bias in the multiplication 
factor of +0.4 Ak/k. Third, for the same level of quadrature order, 
an x-y geometry representation can be used to obtain a much more accurate 
value of the multiplication factor than can be obtained in R-0 geometry. 
Ii2. Recommendations 

Since this work represents the first successful singular eigen- 
function solution of the transport equation for critical cylinders, the 
opportunity exists for much worthwhile additional work in this area. 

For example, a two energy group solution would be very useful for the 
analysis of hydrogen moderated systems. Also, a solution including 
anisotropic scattering could be used to analyze systems in which this 

phenomenon is important. These features have been successfully incorpor- 

29 30 31 

ated into singular eigenfunction solutions ’ ’ of the Boltzmann 

equation in critical slab geometry. How these features might be 
incorporated into the present solution, which is based on the transformed 
integral transport equation, has not been investigated. However, it 
appears to be relatively straightforward to apply the techniques developed 
in the present work to obtain solutions for multiregion problems. The 
results from such solutions could be used as analytic standards for 
methods used to analyze radially fuel- zoned cores or deep radiation 
penetration into multilayered shield configurations. 
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From the results of Chapter VI, it is recommended, that further 
design analysis of the Advanced Power Reactor concept using the DOT-IIW 
discrete-ordinates transport program in the approximation and in R- 9 
geometry include a reactivity bias of -0.4$ Ak/k applied to the calculated 
multiplication factor. It is also recommended that x-y rather than R-0 
geometry be used in analysis requiring high precision in the multi- 
plication factor. 

Finally, the results of Chapters V and VI have shown a serious 
deficiency of low quadrature order discrete-ordinates transport analysis 
in R-0 geometry. Investigation into the source of this problem would 


be a worthwhile area for further research. 
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APPENDIX A 

ORTHOGONALITY AND NORMALIZATION RELATIONS 
FOR THE PSEUDO EIGENFUNCTIONS 


The purpose of this appendix is to present the development of an 
orthogonality relation for the pseudo eigenfunctions found in Chapter III 
and the development of normalization relations for their discrete and 
continuum modes. 

To demonstrate the orthogonality properties, we begin by substituting 
Eq. (3.18) for T] (n>) into the Eq. (3-l6), obtaining 


\ ft*) d(i> = 1. 


(A.l) 


Making the same substitution into Eq. (3.17)> we have 
(v 2 - p 2 )ti v G 0 = cv 2 . 


(A. 2) 


Next we rewrite Eq. (A. 2) for two eigenvalues, v and v 1 , in the form: 


( X “ h)\ = 


(A. 3a) 


('■i-h 


t (M* ) = c. 


(A. 3b) 


Now we operate on Eq. (A. 3a) with the integral of (lO over p and 
operate on Eq. (A. 3b) with the integral of (n) over p, obtaining 


j ^1 - \ (p)\i (M-)dM- = c j T^,(|j,)d|X, 


(A. 4a) 


1 2 1 

J - ^ 2 ) = c J ynjdp.. 


(A.l*) 
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We observe from Eq. (A.l) that the right hand sides of Eqs. (A.4a,b) 
are just equal to c. With this observation, we subtract Eq.(A.4b) 


from Eq. (A. 4a) to obtain 


( ir* - - 7 ) / 


, (m-) = o. 


(A. 5) 


Eq. (A. 5) implies the orthogonality relation 


P 


! \C|i)\ I ^)du = 0, V 4 v. 


(A. 6) 


Next we wish to develop the normalization integral for the 


discrete modes, defined by 


J 1 


p)dp. 


Substituting Eq. (3.19) for (p.) into Eq. (A. 7 ) we have 


f 1 2 cv 2 cv 2 

0 1 " 5-2 ~ 2 ~~ 2 . 

0 V -LX V -LX 

o o 


(A. 7) 


(A. 8) 


By partial fraction expansion Eq. (A. 8) can be shown to be equivalent to 


2 3 r 1 r 

c v o_ f 

+ *^0 


(v Q -u) (v ( 


■ 

d|x . 

\2 

'o + ^ J 


(A. 9) 


Performing the integral over the terms in Eq. (A. 9 ) we obtain 


2 3 

N = o 

° ~ 


[ In I v -M-l 

o 


+ _o__ - in|v +-n| - » 

V -LX 

O V + LL 

0 J o 


which is evaluated at the limits to yield 


N = 

° 2 
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Anticipating the requirement for a discrete normalization integral 
to apply in solving the two-region problem of Chapter IV, we next seek 
the relation 


N 


ol2 


- f 1 

J o 


(n)djj-, 


(A. 12) 


where c > 1, v q1 imaginary; c g < 1, v q2 real. Substituting the discrete 
mode eigenfunctions into Eq. (A. 12) we have 

1 n n 


N 


012 


■/. 


2 c, v 
M- l 1 ol 


I 2 2 

ivl + 


°2 V o2 


V o2 ‘ “ 


djj, # 


(A.15) 


Upon partial fraction expansion Eq. (A. 13) becomes 


N 


ol2 


I (2 2 

1 C 2 | V pl * V o2 f 

V I 2 + v 2 J 

ol 1 o2 c 


o2 


ol 1 


V 62 " ^ 


2 k 

V + ^ 


dp.. 


(A.14) 


Performing the integral over the terms in Eq. (A.l4) , we obtain 


I 1 2 2 

N 10 = C l C 2f V oll V o2 


I 2 2 

l v 0l + v 


N 


ol2 


ol 


V n 1 

V _ + LL 1 

o2 ^ 

2 1 

V - - IX * 
o2 p 

the limits 

to yield 

2 

v o2 fl_ 

i_l • 

2 L c 2 

o2 L 



- ly in 

2i 


Ki I + ^ 


,1 


ol 

V 


- in 


(A.15) 


(A.16) 


Finally, we wish to obtain the normalization function for the 

continuum modes, defined by 
1 

n(v) = f n"H.(n)i\.,(n)dn, v’ = v. 


(v) = I n 2 Ti v (ia.)Ti vl (n)dn, 

J o 


(A.17) 


This function is to be used in the evaluation of coefficients, A(v'), 
appearing in the expansion of an arbitrary function, f(n) in the 
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following manner : 


f(l*) =f ^(v’) ^,(n)dv'. 0 

✓o 

Operating on Eq. (A.l8) by the integral of p T^(|x) over p, we obtain 
1 11 

J M- 2 T| v (M,)f(p- )tln, = j |i 2 Tl v (|x)dp, f A(v* )T|^, (p. )dv‘ (i 

J o J o •'o 

Initially we define the normalization function, N(v), such that 


A(v)N(v) = 


M< )f (m- )d(x 


-f * 2 V 


Wow, substituting Eq. (A.20) into Eq. (A. 19), we obtain 

/ M* A(v»)\,(n)dv’. (A. 21 

o * ° 

Recognizing that the order of integration over singular functions is 
not directly interchangeable, we evaluate Eq. (A.2l) by considering 

1 X 

N(v) = J dv*A(v')J ^ 2 T1 v , (^)Tl v (n)dii. (A. 22 ; 

o 

Upon substitution of Eq, (3*20) for T[^(|Ji) and , (jjl ) , Eq, (A. 22) becomes 
A(v)N(v) J dv'Afv 1 ) rpL^dfij X( v)5(v-m,)\(v ! )&(v'-M* ) 

Jo { 


+ \(v f )&(v , “p,)cP 


2 2 
V - \1 


+ \(v)8(v-|ji)cP v' 


,2 2 
V' - Pi 


+ cP v 


2 2 
v - m- 


,2 2 . 

v 1 - [i 


(A. 23) 
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Evaluating the terms in Eq. (A. 23) which contain delta functions we have 

.1 


A(v)N(v) = A(v)v 2 \ 2 (v) + v 2 cP f dv'A(v') v ,2 \(v') 

Jo 2 ,2 

V - V' 


+ V 


2 2 
+ c v 


r 1 2 

\(v)cP f dv'A(v') y’ 

J ,2 

o v' - 

J" dv'A^v' )v ,2 ^'V 


V' - V 
2 


d^i , 


(A. 24) 


2 ,2 2 

- U v' - (J- 


The last term in Eq. (A. 24) is expanded by partial fractions to yield 

1 

c 2 v 2 / dv'A(v' ) v'‘ 

Jr> .? 


,2 2 

v ! - V 


t\ 


P V 


P v* 


dj-L . 


2 2 ,2 2 

v ~ M- v' - p, 


Next we use the function, \(v), in the form obtained by integrating Eq. 
(3.20) over n, 


\(v) = 1 - cP 
to write this term as 


j 


$ 2 
o v - p 


<4v, 


cv 


f dv'A(v' ) V * 2 [\(v») - \(v)] . 

J ,2 2 

~ v' - v 


The term as written above identically cancels the second and third 
terras of Eq. (A. 24), leaving 
A(v)N(v) = A(v)v 2 X 2 (v). 

Now returning to the evaluation of Eq. (A.2l) we insert the 
eigenfunctions to obtain 


(A.25) 


A(v)N(v) = 


v) =J 1 du^ 2 | 


cP v + \(v)8 (v-ji )[ I dv'A(v' He P 


2 2 
v 


ir- 


,2 2 

v’ — |J, 


+ \(v' )6(v'-n) | • (A. 26) 
After performing the integrals over the terms containing delta functions 
in Eq. (A. 26) we have 
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A(v)N( v ) = A(v)v 2 \ 2 (v) + cv 2 f dv'A(v')P v' g \(v' ) 
1 Jo 2 

\(v) I dv'A(v' ) P 

Jo .2 


2 ,2 

v - v' 


+ cv 


,2 2 

V 1 -V 


cP v 2 p 2 dp f dv'A(v')cP v ' 2 (A. 27) 

2 2 •'o .22 

v -p v' -p 

We note that the first three terms of Eq. (A. 27) are identical to those 

of Eq. (A. 24). The fourth terms differ by the order of integration. 

We evaluate the fourth term of Eq. (A. 27 ) by putting it in the form of 

28 

the Poincare' -Bert rand formula as follows: 



c 2 P 


>f u 2 v g dppf A(y') v 1 

J 0 2 2 Jo .2 


dv' 


2 

= -c v P 


V -p 
1 


,2 2 

v' -p 


pf .d^ PjJ g(v',p) dv* 

I v* - p 

Jo p-v 

2 2 

where g(v',p) = A(v' )v' p /(v+p)(v'+p) . 
The Poincare' -Bertrand formula is 


(A. 28) 


p Jo g( v » dv’ 
p-v v'-p 

= -TT 2 g(v,v) + dv' P 1_ P g(v',p)dp. (A. 29 ) 

V-ji V 1 -jX 

Now, comparison of Eqs. (A. 28), (A. 29 ) and the fourth term of Eq. (A. 24) 

shows that the difference of the two fourth terms is the additional term, 
2 2 2 

c v TT g ( v, v) , arising from Eq. (A. 29). With this observation, we cancel 
the identical terms of Eq. (A. 27 ) as before and write Eq.(A.27) as 

A(v)N(v) = A(v)v 2 A 2 (v) + c 2 v 2 rr 2 g(v,v). (A. 50 ) 

Finally, substituting Eq. (A. 28) for g(v,v) into Eq. (A. 30 ) and canceling 


like terms we obtain 



as the normalization function for the continuum modes of the pseudo 


eigenfunctions 
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APPENDIX B 

PSEUDO EIGENFUNCTION COMPLETENESS 

Theorem: The functions T^(p) and T^(u-) form a complete set in the sense 

that any arbitrary known even function $(p) can be expanded in terms of 
these eigenfunctions for -1 < u < 1. Proof: We write the expansion 

*(n) = + ti (B.l) 

and the task is to show that we can solve for the expansion coefficients 
a^, A(v), in terms of the known function $(p). 

In the proof we make use of the boundary values of the dispersion 
function,, A(z), where 

A - (v) = 2 + cv 2 Pf 1 , du ± nicv, 

J “ 1 2 2 
LL “ V 

or r 

A” (v) = \(v) ±TTiCV, 

28 

are obtained from the Plemelj formulas 

F* (y) = Pj^ f(x) dx ± nif(y). 
x-y 

We note that A(v) and \(v) are twice as large as the functions defined 
in Eqs . (2.13) and (2.15). This difference is consistent with the 
half range normalization to unity, Eq. (3.16). First we write from 
Eq. (B.l) 

#'G0 = A ( v ) \ (^) d v (B.2a) 

where 

*'(n) = *(m0 - 

Inserting Eq. (3*20 ) for T| v (m-) into Eq. (B.2a) we have 


(B.2b) 
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# ' GO = XGQAGO + c P J 1 A(v)v 2 dv, 

2 2 2 2 
V - [JL 


(B.3a) 


2 

*'(h)=t (a + (p) + a‘M)a(iO +f Pj\ A | v)v 3 dv, (B.3b) 


or 

2 
y) 

r 

V — jir^ 

where we have substituted the sum of the boundary values of the 
dispersion function for 

Next we introduce the function N(z) defined by 


N(z) = W Jti 


cv A(v) 

2 , 2 
V - z 


dv . 


(B.4) 


We note that N(z) e A in the complex plane cut from -1 to 1 along the 
real axis and N(z) ~ l/z as z-"=°. Adding and subtracting the boundary 
values of N(z) as given by the Plemelj formulas, we obtain 

c P r 1 , A(v)v 2 dv = ttl(n + (Q + H*(n) ), (B.5a) 

o 0-1 o o \ ' 


2 2 
v - M- 


and 


AGO = ~ (n + (h) - if GO )• 


(B.5b) 


Substitution of Eqs . (B.5a,b) into Eq. (B.3b) yields 


^ * Cm-) = ^ {(a + (.) + a*(m)) (n + {^) - n”G0) 

+ 2rric|i. (|i, ) + N (|0^ j • (B .6 ) 

Now, inserting the difference of the dispersion function boundary 
values for 2Tric|j,, multiplying through and canceling like terms, Eq. 
(B.6) becomes 

cm-* * (m) = a + G0n + G0 - a’ (On" GO* (b . 7 ) 

Next we seek a function, N(z)A(z), whose boundary values on the cut 
satisfy equation (B.7)» Thus we define the function <J(z) as 
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J(z) = N(z)A(z) - c J 1 . 'P 2 § ' Cm. ) dp, 

r>, * ” -L o n 


(B.8) 


2ni 


2 2 
p - z 


with boundary values such that 

J + (p) - J”(|i) = N + (p)A + (p) - flf(p)A"(p) - cp«*(p). (B.9) 

Comparison of Eq. (B.9) and Eq. (B.7) shows that the difference of the 
J(z) function boundary values is zero. With no discontinuity across 

the cut, J(z) is a bounded entire function which, by Liouville's 

2 

theorem, must be a constant. Furthermore, since j(z)^o as z J(z) 
is everywhere equal to zero. With this result we rewrite Eq. (b. 8) as 

2 

1 t 


H(z) = 2rtA(z) £l d “» hi 2 


(B.10) 

p“ - z~ 

Now since ± v q are the roots of the dispersion function, in order for 
N(z) as defined in Eq. (B.10) to be analytic at ± v , we must require 


%Ll l in . ) 

x o o 


2 2 ^ = °‘ 


(B.lla) 


To obtain the solution for the discrete coefficient a , we insert 

o' 

Eq. (B.2b) into Eq. (B.lla), yielding 


1-1 -Sjy 5 ^ = °* 

p - Vq 


Substituting the definition of T^(p), Eq. 
the full range of v and p, we rewrite Eq. 


(3.19)> which 
(B.llb) as 


(B.llb) 
is valid for 


-1 M- 2 T^(p)^(p) - a Q 1^(10,)) dp = 0. (B.llc) 

v 2 

o 

Now use of the normalization relation Eq. (A. 7) yields the desired 
expression for a Q , namely 



Ill 


a 0 = S' ti 


(B.12) 


To obtain the similar expression for the continuum coefficient 
A(v), we insert the boundary values of N(v) as defined from Eq. (B.IO) 
into (B.5b) to give 


•i n-S y) a(1 , 


ACu) = -2_i^_r^ ^ 

' 2ni |_A + U0 " A’(n)J _1 H' 2 - 


+ 1 


sk iLM 

A + (p,) A“(n) I 2 


(B.13a) 


Upon substitution of the sum and difference of the dispersion function 
boundary values Eq. (B.13a) becomes 

A(t*) = 2 {p jt-L ^ ,2 §‘ (n*) dp,* + i . (B.13b) 


A (p-)A (p) 


M- - M- 


Now letting \i = v and p, 1 = jj, we have 


A(v) = 


v 2 A + (v)A (v) 


P £ cvV(^ 2 dp, + £ v 2 \(v)$ 1 (v)8(v-p,)dp,\ . 

m 2 _ „2 " j 


v" - p, 


(B.13c) 

The two terms of the continuum eigenfunction as defined by Eq. (3 ,20) 
appear in Eq. (B.13c), allowing us to write 


A(v) = 2 


v 2 A + ( v ) A** ( v ) 


J^l M' 2 T]^(p,)$ 1 (n)dp,. 


(B.13d) 


Now inserting Eq. (B.2b) for f(p,) into Eq. (B.13d) we obtain 

A(v) = -5^7 : — 

v A (v)A (v) 


jij_ ^^(p.) (§ ’ Cm-) - a 0 T l 0 (P‘) ) d ^- 


(B.13e) 


Applying the orthogonality relation, Eq. (A.6), and substituting the 
normalization constant N(v), Eq. (A.Jl), we have the solution for A(v) as 
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A(v) = h\GO*0O^. (B.14) 

This completes the proof for the full range case. 

An immediate extension can be made to the half range case if we 
note from Eqs. (3.17) and (3.18) that 

and from the theorem, $(-|0 = 5 GO* 

Thus we can write the solution for A(v) as 

Mv)=5TO To (B ' 15) 

and the expansion given by Eq. (B.l) as 

*60 = A(v)\(Odv, (B.16) 

where, again from Eq. (3*17) and (3*l8), we have used 

T) (n) = T)^^) and we have assumed, because of the evenness of 
all functions involved, that A(-v) = A(v). 
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APPENDIX C 

NEUTRON DENSITY DISTRIBUTIONS FOR REFLECTED CYLINDERS 

Neutron density distributions corresponding to several of the cases 
for which the critical dimensions are given in Table V are presented 
here for the use of interested researchers. Since the numerical results 
of all one hundred cases calculated would be rather voluminous, we 
have selected combinations of cases which provide various parametric 
ranges. First, we present the largest and the smallest cases calculated. 
Then follows a variety of cases from which the effects of the independent 
variation of core c value, reflector c value and reflector thickness 


can be drawn 
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c^=l . 02 , 

c 2 =0.99 

c^l. 4 , c 

2=0.85 

or 

R 2 " R 1 = 

20 MFP 

R 2 " R l = 

1 MFP 

f 2 =r - R l 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

V R i 





0.0 

1.0000 

0.5415 

1.0000 

0.4906 

0.01 

0.9999 

0.5126 

0.9999 

0.4788 

0.02 

0.9998 

0.4871 

O.9998 

0.4691 

0.1 

0.9949 

0.3327 

0.9946 

0.4082 

0.2 

0.9795 

0.2124 

0.9785 

0.3506 

o.3 

0.9542 

O.I378 

0.9520 

0.3037 

0.4. 

0.9192 

0.0902 

0.9153 

0.2640 

0.5 

0.8752 

0.0593 

0.8688 

0.2294 

0.6 

0.8228 

0.0387 

O.8131 

0.1987 

0.7 

0.7628 

0.0247 

0.7485 

0.1710 

0.8 

0.6959 

0.0149 

0.6753 

0.1455 

0.9 

0.6228 

0.0076 

0.5927 

0.1212 

0.95 

0.5838 

0.0046 

0.5464 

0.1089 

O.98 

0.5592 

0.0029 

0.5155 

0.1010 

0.99 

0.5507 

0.0022 

0.504l 

0.0981 

1.00 j 

0.5415 

0.0015 

0.4906 

0.0948 
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r x =r /^i 

or 

c^=l .02, 
R 2 ' R 1 = 

c 2 = 0.95 
3 MFP 

c =1.02, 
R 2 - R 1 = 

= 2 = 0.95 
10 MFP 

r 2 =r " K l 

p(*\ )/p(o) 

p(r 2 )/p(o) 

p(r x )/ p(o) 

p(r 2 )/p(o) 

R -R_ 

-L 




£ 1 





0.0 

1.0000 

O.2858 

1.0000 

0.3153 

0.01 

0.9999 

0.2806 

0.9999 

0.2990 

0.02 

0.9996 

O.2758 

0.9997 

0.2847 

0.1 

0.9912 

0.2423 

0.9916 

0.1977 

0.2 

0.9649 

0.2068 

0.9668 

0.1277 

0.3 

0.9219 

0.1761 

0.9261 

0.0830 

0.4 

0.8633 

0.1490 

0.8705 

0.0542 

0.5 

0.7906 

0.1249 

0.8oi4 

0.0354 

0.6 

0.7058 

0.1031 

0.7206 

0.0230 

0.7 

0.6110 

0.0834 

0.6300 

0.0147 

0.8 

0.5085 

0.0651 

0.5316 

0.0090 

0.9 

o.4oo6 

0.0476 

0.4274 

0.0050 

0.95 

0.3450 

O.O387 

0.3735 

0.0033 

0.98 

0.3106 

0.0330 

0.3397 

0.0024 

0.99 

0.2987 

0.0309 

0.3280 

0.0020 

1.0 

0.2858 

0.0283 

0.3153 

0.0016 










3 


0^=1. 05, c 2 =0.95 


R 2 “ R i = 5 MFP 


c x =1.05, c 2 =0.95 


R 2 - R x = 10 MFP 


r i =r/E i 

or 

r 2 = r “ R l 

V R 1 

0.0 

0.01 

0.02 

0.1 

0.2 

0.3 

0.4 

0.5 

o.6 

0.7 

0.8 

0.9 

0.95 

0.98 

0.99 


p(r 1 )/p(o) 


1.0000 

0.9999 

0.9997 

0.9931 

0.9726 

0.9388 

0.8926 

0.8547 

0.7665 

0.6887 

0.6052 

0.5108 

0.4617 

0.4507 

0.4197 

0.4076 


p(r 2 )/p(o) 


0.4076 

0.5989 

0.3912 

0.3386 

0.2848 

0.2394 

0.2003 

0.1662 

0.1361 

0.1091 

0.0845 

0.0613 

0.0497 

0.0423 


p(r 1 )/p(o) 


1.0000 

0.9999 

0.9997 

0.9936 

0.9745 

0.9431 

O.8999 

0.8458 

O.7816 

0.7086 

O.6278 

0.5399 

0.4928 

0.4628 

0.4522 


p(r 2 )/p(o) 


o.44o4 

o.4l4i 

0.3918 

0.2622 

0.1641 

0.1042 

0.0668 

0.0429 

0.0275 

0.0174 

0.0106 

0.0058 

0.0039 

0.0027 

0.0023 


1.0 


0.0395 

0.0362 


0.44o4 


0.0019 
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r r r / R i 

^=1.05, 

c 2 =0.85 

Cj=1.05, 

c 2 =0. 8 5 

or 

R 2 - R 1 = 

3 MFP 

E 2 - R 1 = 

10 MFP 

r 2 =r-E 1 

V R 1 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

0.0 

1.0000 

0.2859 

1.0000 

O.2897 

0.01 

0.9999 

O.2765 

0.9999 

0.2624 

0.02 

0.9997 

O.2685 

0.9997 

0.2405 

0.1 

0.9913 

0.2177 

0.991^ 

0.1291 

0.2 

0.9655 

0,1708 

0.9658 

0.0625 

0.5 

0.9233 

0.1349 

0.9239 

0.0309 

0.4 

0.8657 

0.1067 

0.8667 

0.0155 

0.5 

0.7942 

0.0842 

0.7956 

0.0078 

0.6 

0.7105 

0.0660 

0.7125 

o.oo4o 

0.7 

0.6167 

0.0510 

0.6192 

0.0020 

0.8 

0.5148 

0.0384 

0.5179 

0.0010 

0.9 

o.4o64 

0.0275 

0.4100 

0.0005 

0.95 

0.5494 

0.0223 

0.3531 

0.0003 

0.98 

0.3133 

0.0190 

0.3171 

0.00021 

0.99 

0.3004 

O.OI78 

0.3042 

0.00018 


0.2859 


L 


0.2897 


1.0 


0.0164 


0.00015 
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r r r/R i 

C-J=l.l, c 

2 =0.95 

c^=l.l, c 

2=0.95 

or 

R 2 " R 1 = 

3 MFP 

R 2 ' R 1 = 

10 MFP 

r 2 =r-R l 

V R 1 

p(r 1 )/p(o) 

p(r 2 )/p(°) 

p(r 1 )/p(o) 

p (r 2 )/p(o) 

0.0 

1.0000 

0.4992 

1.0000 

0.5295 

0.01 

0.9999 

0.4864 

0.9999 

0.4916 

0.02 

0.9998 

0.4753 

0.9998 

0.4611 

0.1 

0.99^5 

0.4031 

0.9949 

0.2953 

0.2 

0.9780 

0.3329 

0.9796 

0.1789 

0.3 

0.9508 

0.2758 

0.9545 

0.1113 

0.4 

0.9134 

0.2280 

0.9197 

0.0702 

0.5 

0.8662 

0.1871 

0.8758 

0.0446 

0.6 

0.8099 

0.1518 

0.8233 

0.0284 


0.7 

0.8 

0.9 

0.95 

0.98 

0.99 


0 . 7^52 
0.6728 
0.5927 
0.5491 
0.5208 
0 . 5106 
0.4992 


0.1206 

0.7628 

0.0178 

0.0927 

0.6947 

0.0108 

0.0669 

0.6189 

0.0059 

0.0541 

0.5773 

0.0039 

0.0459 

0.5502 

0.0028 

0.0429 

0. 54o4 

0.0024 


1.0 


0.0392 


0.5293 


0.0019 
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r r r/E i 

or 

c =1.2, c, 
± < 

R 2 “ R 1 = 

5 =0.95 
3 MFP 

c^=1.2, c 
R 2 ■ R 1 = 

2 =0.85 
3 MFP 

r a =r ' E i 

V R 1 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

0.0 

1.0000 

0.5732 

1.0000 

0.4639 

0.01 

0.9999 

0.5539 

0.9999 

0.4266 

0.02 

0.9998 

0.5380 

0.9998 

0.3976 

0.1 

0.9956 

0.4413 

0.9942 

0.2483 

0.2 

0.9823 

0-3546 

0.9768 

0.1487 

0.3 

0.9602 

0.2879 

0.9482 

0.0922 

0.4 

0.9297 

0.2342 

O.9087 

0.0582 

0.5 

0.8910 

0.1898 

O.8589 

0.0371 

0.6 

0.8445 

0.1523 

0.7995 

0.0237 

0.7 

0.7903 

O.II99 

0.7310 

0.0151 

0.8 

0.7286 

0.0915 

0.6540 

0.0094 

0.9 

0.6587 

0.0655 

0.5680 

0.0055 

0.95 

0.6196 

0.0528 

0.5204 

o.oo4o 

0.98 

0.5937 

0.0447 

0.4889 

0.0031 

0.99 

0.5842 

o.o4i7 

0.4774 

0.0027 

1.0 

0.5732 

0.0382 

0.4639 

0.0024 
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r^r/Ri 

or 

V r - E i 

V\ 

0 ^= 1 . 4, c 
E 2 ‘ \ = 

2 =o.95 
3 MFP 

c^l.4, c 
R 2 “ R 1 = 

2=0.95 
10 MFP 

p(r 1 )/p(o) 

p(r 2 )/p(o) 

O 

Q_ 

1 — 1 
CL 

p(r 2 )/p(o) 

0.0 

1.0000 

0.4992 

1.0000 

0.6384 

0.01 

0.9999 

0.4864 

0.99996 

0-5597 

0.02 

0.9998 

0.4755 

0.99985 

0.5057 

0.1 

0.9945 

0.4031 

0.9965 

0.2790 

0.2 

0.9780 

0.3329 

0.9859 

0.1564 

0.3 

0.9508 

0.2758 

0.9683 

0.0933 

0.4 

0.9134 

0.2280 

0.9438 

0.0573 

0.5 

0.8662 

0,1871 

0.9128 

0.0357 

0.6 

0.8099 

0.1518 

0.8744 

0.0224 

0-7 

0.7452 

0.1206 

0.8275 

0.0139 

0.8 

0.6728 

0.0927 

0.7774 

0.0084 

0.9 

0.5927 

0.0669 

0.7168 

0.0045 

0.95 

0.5490 

0.0541 

0.6819 

0.0030 

0.98 

0.5208 

0.0459 

0.6581 

0.0021 

0.99 

0.5106 

0.0429 

0.6491 

0.0018 


0.4992 


0.6384 


1.0 


0.0392 


0.0014 














